{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "62845e91",
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.optim as optim\n",
    "import torch.utils\n",
    "import torch.utils.data\n",
    "from torch.utils.data import DataLoader\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from math import pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "d9d302df",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exact solution\n",
    "def exact_solution(x, t):\n",
    "    return torch.sin(x)*torch.cos(4*pi*t)\n",
    "\n",
    "# Initial condition\n",
    "def initial_condition(x):\n",
    "    return torch.sin(x)\n",
    "\n",
    "def initial_condition_time(x):\n",
    "    return 0*np.pi**2*torch.sin(np.pi*x)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "cbd4c671",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<torch._C.Generator at 0x7f0ba400c790>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# assigning number of points\n",
    "initial_pts = 2000\n",
    "left_boundary_pts = 2000\n",
    "right_boundary_pts = 2000\n",
    "residual_pts = 10000\n",
    "\n",
    "# Type of optimizer (ADAM or LBFGS)\n",
    "opt_type = \"LBFGS\"\n",
    "manualSeed = 1\n",
    "\n",
    "#np.random.seed(manualSeed)\n",
    "#random.seed(manualSeed)\n",
    "torch.manual_seed(manualSeed)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "61ba54bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# initial points\n",
    "x_init = pi*torch.rand((initial_pts,1)) # initial pts\n",
    "t_init =  0*x_init\n",
    "init =  torch.cat([x_init, t_init],1)\n",
    "u_init = initial_condition(init[:,0]).reshape(-1, 1)\n",
    "u_init_t = initial_condition_time(init[:,0]).reshape(-1, 1) #new"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "1334fc49",
   "metadata": {},
   "outputs": [],
   "source": [
    "#boundary points\n",
    "\n",
    "xb_left = torch.zeros((left_boundary_pts, 1)) # left spatial boundary\n",
    "tb_left = torch.rand((left_boundary_pts, 1)) # \n",
    "b_left = torch.cat([xb_left, tb_left ],1)\n",
    "u_b_l = 0*xb_left\n",
    "u_b_l_xx = 0*xb_left #new\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "fc482504",
   "metadata": {},
   "outputs": [],
   "source": [
    "xb_right = pi*torch.ones((right_boundary_pts, 1)) # right spatial boundary\n",
    "tb_right = torch.rand((right_boundary_pts, 1)) # right boundary pts\n",
    "b_right = torch.cat([xb_right, tb_right ],1)\n",
    "u_b_r = 0*xb_right\n",
    "u_b_r_xx = 0*xb_right #new\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "1d80068f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# collocation/ interior points\n",
    "x_interior = pi*torch.rand((residual_pts, 1))\n",
    "t_interior = torch.rand((residual_pts, 1))\n",
    "interior = torch.cat([x_interior, t_interior],1)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "d2a71073",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Training set\n",
    "training_set = DataLoader(torch.utils.data.TensorDataset(init, u_init, u_init_t, b_left, b_right), batch_size=2000, shuffle=False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "415db5b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# neural network Class\n",
    "\n",
    "class NeuralNet(nn.Module):\n",
    "\n",
    "    def __init__(self, input_dimension, output_dimension, n_hidden_layers, neurons):\n",
    "        super(NeuralNet, self).__init__()\n",
    "        # Number of input dimensions n\n",
    "        self.input_dimension = input_dimension\n",
    "        # Number of output dimensions m\n",
    "        self.output_dimension = output_dimension\n",
    "        # Number of neurons per layer \n",
    "        self.neurons = neurons\n",
    "        # Number of hidden layers \n",
    "        self.n_hidden_layers = n_hidden_layers\n",
    "        # Activation function \n",
    "        self.activation = nn.Tanh()\n",
    "        \n",
    "        self.input_layer = nn.Linear(self.input_dimension, self.neurons)\n",
    "        self.hidden_layers = nn.ModuleList([nn.Linear(self.neurons, self.neurons) for _ in range(n_hidden_layers)])\n",
    "        self.output_layer = nn.Linear(self.neurons, self.output_dimension)\n",
    "\n",
    "    def forward(self, x):\n",
    "        # The forward function performs the set of affine and non-linear transformations defining the network \n",
    "        # (see equation above)\n",
    "        x = self.activation(self.input_layer(x))\n",
    "        for k, l in enumerate(self.hidden_layers):\n",
    "            x = self.activation(l(x))\n",
    "        return self.output_layer(x)\n",
    "# Model definition\n",
    "my_network = NeuralNet(input_dimension = init.shape[1], output_dimension = u_init.shape[1], n_hidden_layers=4, neurons=20)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "da5ea1e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def init_xavier(model, retrain_seed):\n",
    "    torch.manual_seed(retrain_seed)\n",
    "    def init_weights(m):\n",
    "        if type(m) == nn.Linear and m.weight.requires_grad and m.bias.requires_grad:\n",
    "            g = nn.init.calculate_gain('tanh')\n",
    "            torch.nn.init.xavier_uniform_(m.weight, gain=g)\n",
    "            #torch.nn.init.xavier_normal_(m.weight, gain=g)\n",
    "            m.bias.data.fill_(0)\n",
    "    model.apply(init_weights)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "6bf0ace9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Random Seed for weight initialization\n",
    "retrain = 128\n",
    "# Xavier weight initialization\n",
    "init_xavier(my_network, retrain)\n",
    "\n",
    "if opt_type == \"ADAM\":\n",
    "    optimizer_ = optim.Adam(my_network.parameters(), lr=0.001)\n",
    "elif opt_type == \"LBFGS\":\n",
    "    optimizer_ = optim.LBFGS(my_network.parameters(), lr=0.1, max_iter=1, max_eval=50000, tolerance_change=1.0 * np.finfo(float).eps)\n",
    "else:\n",
    "    raise ValueError(\"Optimizer not recognized\")\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "d8ce5a43",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fit(model, training_set, interior, num_epochs, optimizer, p, verbose=True):\n",
    "    history = list()\n",
    "    \n",
    "    # Loop over epochs\n",
    "    for epoch in range(num_epochs):\n",
    "        if verbose: print(\"################################ \", epoch, \" ################################\")\n",
    "\n",
    "        running_loss = list([0])\n",
    "        \n",
    "        # Loop over batches\n",
    "        for j, (initial, u_initial, u_initial_t, bd_left, bd_right) in enumerate(training_set):\n",
    "            \n",
    "            def closure():\n",
    "                # zero the parameter gradients\n",
    "                optimizer.zero_grad()\n",
    "                \n",
    "                initial.requires_grad = True\n",
    "                 # for initial\n",
    "                u_initial_pred_ = model(initial)\n",
    "                inputs1 = torch.ones(initial_pts, 1 )\n",
    "                grad_u_initial = torch.autograd.grad(u_initial_pred_, initial, grad_outputs=inputs1, create_graph=True)[0]\n",
    "                u_initial_pred_t_ =  grad_u_initial[:, 1]\n",
    "                \n",
    "                \n",
    "                # boundary\n",
    "                bd_left.requires_grad = True\n",
    "                bd_right.requires_grad = True\n",
    "                u_bd_left_pred_ = model(bd_left)\n",
    "                u_bd_right_pred_ = model(bd_right)\n",
    "                inputs2 = torch.ones(left_boundary_pts, 1)\n",
    "                inputs3 = torch.ones(right_boundary_pts, 1)\n",
    "                grad_u_b_l = torch.autograd.grad(u_bd_left_pred_, bd_left, grad_outputs=inputs2, create_graph=True)[0]\n",
    "                grad_u_b_r = torch.autograd.grad(u_bd_right_pred_, bd_right, grad_outputs=inputs3, create_graph=True)[0]\n",
    "                u_b_l_x = grad_u_b_l[:, 0]\n",
    "                u_b_r_x = grad_u_b_r[:, 0]\n",
    "                u_b_l_xx = torch.autograd.grad(u_b_l_x, bd_left, grad_outputs=torch.ones(bd_left.shape[0]), create_graph=True)[0]\n",
    "                u_bd_left_pred_xx_ = u_b_l_xx[:, 0]\n",
    "                \n",
    "                u_b_r_xx = torch.autograd.grad(u_b_r_x, bd_right, grad_outputs=torch.ones(bd_right.shape[0]), create_graph=True)[0]\n",
    "                u_bd_right_pred_xx_ = u_b_r_xx[:, 0]\n",
    "                \n",
    "                \n",
    "                \n",
    "                \n",
    "                \n",
    "                \n",
    "               \n",
    "                # residual calculation\n",
    "                interior.requires_grad = True\n",
    "                u_hat = model(interior)\n",
    "                inputs = torch.ones(residual_pts, 1 )\n",
    "              \n",
    "                grad_u_hat = torch.autograd.grad(u_hat, interior, grad_outputs=inputs, create_graph=True)[0]\n",
    "                \n",
    "                u_x = grad_u_hat[:, 0]\n",
    "                u_t =  grad_u_hat[:, 1]\n",
    "                \n",
    "                grad_grad_u_t = torch.autograd.grad(u_t, interior, grad_outputs=torch.ones(interior.shape[0]), create_graph=True)[0]\n",
    "                u_tt = grad_grad_u_t[:, 1]\n",
    "                \n",
    "                grad_grad_u_x = torch.autograd.grad(u_x, interior, grad_outputs=torch.ones(interior.shape[0]), create_graph=True)[0]\n",
    "                u_xx = grad_grad_u_x[:, 0]\n",
    "                \n",
    "                grad_grad_u_xx = torch.autograd.grad(u_xx, interior, grad_outputs=torch.ones(interior.shape[0]), create_graph=True)[0]\n",
    "                u_xxx = grad_grad_u_xx[:, 0]\n",
    "                \n",
    "                grad_grad_u_xxx = torch.autograd.grad(u_xxx, interior, grad_outputs=torch.ones(interior.shape[0]), create_graph=True)[0]\n",
    "                u_xxxx = grad_grad_u_xxx[:, 0]\n",
    "                \n",
    "            \n",
    "                \n",
    "                # Item 1. below\n",
    "                loss_ini = torch.mean((u_initial_pred_.reshape(-1, ) - u_initial.reshape(-1, ))**p) + torch.mean((u_initial_pred_t_.reshape(-1, ) - u_initial_t.reshape(-1, ))**p)\n",
    "                \n",
    "                loss_p = torch.mean((u_tt.reshape(-1, ) + u_xxxx.reshape(-1, )-(1 -16*pi**2)*torch.sin(interior[:, 0])*torch.cos(4*pi*interior[:, 1]))**p)\n",
    "                \n",
    "                loss_bd = torch.mean((u_bd_left_pred_.reshape(-1,))**p) + torch.mean((u_bd_right_pred_.reshape(-1,))**p)+torch.mean((u_bd_left_pred_xx_.reshape(-1,))**p) + torch.mean((u_bd_right_pred_xx_.reshape(-1,))**p)\n",
    "                \n",
    "                loss = loss_ini + 0.1*loss_p + loss_bd\n",
    "               \n",
    "                # Item 2. below\n",
    "                loss.backward()\n",
    "                # Compute average training loss over batches for the current epoch\n",
    "                running_loss[0] += loss.item()\n",
    "                return loss\n",
    "            \n",
    "            # Item 3. below\n",
    "            optimizer.step(closure=closure)\n",
    "            \n",
    "        print('Loss: ', (running_loss[0] / len(training_set)))\n",
    "        history.append(running_loss[0])\n",
    "\n",
    "    return history\n",
    "n_epochs = 15000\n",
    "#history = fit(my_network, training_set, interior, n_epochs, optimizer_, p=2, verbose=True )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "8cfaebdc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "NeuralNet(\n",
       "  (activation): Tanh()\n",
       "  (input_layer): Linear(in_features=2, out_features=20, bias=True)\n",
       "  (hidden_layers): ModuleList(\n",
       "    (0): Linear(in_features=20, out_features=20, bias=True)\n",
       "    (1): Linear(in_features=20, out_features=20, bias=True)\n",
       "    (2): Linear(in_features=20, out_features=20, bias=True)\n",
       "    (3): Linear(in_features=20, out_features=20, bias=True)\n",
       "  )\n",
       "  (output_layer): Linear(in_features=20, out_features=1, bias=True)\n",
       ")"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# saving and loading Model\n",
    "FILE = \"first_.pth\"\n",
    "#torch.save(my_network, FILE)\n",
    "\n",
    "# uncomment below when you need to test for different points\n",
    "my_network = torch.load(FILE)\n",
    "my_network.eval()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "7f29e6b2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Relative Error Test:  0.0005351814706955338 %\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAEYCAYAAACKvFuOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABEd0lEQVR4nO2de3weVZ3/39/e0wtNL7SlLW1KKdTlfpObiKW7BazuFhVcWaqAW1RgEVd+Sy24iAu0XlAQcRUR64K39fdDUCoCikWgApbYEhtC00CaNo0NpWlD2jRN8nx/f8w8IUmf28wz88w56Xm/Xs9rkrmc+Xye7zzznTNz5hxRVRwOh8PhiJNBSQtwOBwOx8DHJRuHw+FwxI5LNg6Hw+GIHZdsHA6HwxE7Ltk4HA6HI3aGJC3ARCZOnKgVFRVJy3A4HA6rePnll3eo6qGZlrlkk4GKigrWrl0batu6ujpmz54dsaLSYrsH2/WD/R5s1w/2e0hCv4hszrbM3UaLmPHjxyctoWhs92C7frDfg+36wX4Ppul3ySZi9u7dm7SEorHdg+36wX4PtusH+z2Ypt8lm4gZNMj+r9R2D7brB/s92K4f7Pdgmn4jntmIyEeAc4ETgROAMcCPVfWyEGVNB74MXABMAJqAR4BbVbUlIslZGTp0aNy7iB3bPdiuH97xULF0VUHrr/z4ibzv76bFKSkQAykGtmKafiOSDXAzXpJpA7YCc8MUIiKzgTXAJOBRoAZ4N/BZ4AIROVtV3ypGaGtrK83NzXR2dmZc3tXVxZtvvlnMLhLHZg9Dhw4llUoxceLEpKUEon9SOXNSij81F35levn/rAPWHTD/3n8+joUnzihOXAja2tqsi0F/bPdgmn5Tks3n8JLMJrwazh9ClvMdvERznarek54pIt/w93E78OmwIltbW9m+fTvTpk2jrKwMETlgne7ubgYPHhx2F0ZgqwdVpb29nS1bttDa2sohhxyStKSs5KuxbGo98NgKwzU/q+Kan1X1/F+/YmEk5ebDpJNcWGz3YJp+I27qqeofVLVWi+iC2q/VLADqgXv7Lb4F2AMsFpFRYffR3NzMtGnTGDlyZMZEA7B///6wxRuDrR5EhJEjRzJp0iSam5uTlnMAFUtX9XzyccrEeHpjD6KhGLZu3Rpr+aXAdg+m6TelZhMF8/zpk6qa6r1AVd8WkefxktEZwO/D7KCzs5OysrKc64wYMSJM0UZhu4exY8fyt7/9LWkZQOHPXPrz9LZoaja56K0t6hrPkUceGWl5SWC7B9P0G1GziYij/enGLMtr/elRmRaKyFUislZE1jY1NbFjxw6amppobGykpaWFuro6VJV9+/ahquzZswegz1RVaWtrI5VKsW/fPrq6uti/fz/79++ns7OTjo4Ouru7aW9vR1V7mib2L2vv3r09ZXR3d9PR0UFnZ2efMvbt20cqlcpZRvq2Uu8y0nq6urp6ykjrSW/79ttv9/HU3t5ulaf29nY6OztJpVJUVlYC8PLLLwNQWVlJKpWiurqa9vZ26urqaGlpobGxkXTc6+vraWtro6amhq6uLtavX9+njPS0qqqKjo4OamtraW1tpaGhgebmZv7l7se5+M5fc+aXH+O8qSlGDVEumtkNwOIj+04vntXNsEHKBdNTTByunDA+xaQRytxypXyYMnyQUjZYEZThgxVQhg7yaj1jhnrTKWXe/LljU4wcopwyIcXho5STJqQ4dlyK2WOUMyelOHSEt59hg5SLZ/XV8cX7f0XF0lV87r9/ldFTc3MzDQ0NtLa2UltbS0dHB1VVVRm/l/Xr11NVVUVNTQ1tbW3U19dn/D21t7dTXV2dWJzyedqwYUMfT11dXVZ5WrNmTd44Re0pF2La4Gki8j68ZzaBWqOJyH3AEmCJqt6fYfntwDJgmaouz1XWqaeeqpl6EHj11Vd517veVagkR4IkEau4b02VmlI933EMHETkZVU9NdOygVSzMYL0VbjNZPJw+eWXc/PNNyegJjiljkEcz0D+ZXZ3pOWFoRhf6Stfm7Hdg2n6B1Ky2e1Px2ZZnp6/K04Ro0aFbn9QNBUVFZSVlTF69GgmT57M5ZdfTltbW+ByivVQTGL69re/zamnnsrw4cO5/PLLQ5VRqhjE+aD9x3XmtAYM4/OUU06JSU3psN2DafoHUrJ5zZ9mfCYDzPGn2Z7pRELSNZtf//rXtLW1UVlZydq1a7ntttsOWKerqytnGUl6mDp1KjfffDNXXnll6DLi1v/he/8Y+y2z9LMUkwiSdEy7qg6D7R5M0z+Qkk363ZwFItLHl4iMAc4G9gIvxCkiyZpNb6ZNm8aFF17IX//6V8BrFnzvvfcyZ84c5szx8u5jjz3GiSeeSHl5OWeddRavvPIK4Hn4y1/+wsknn8yYMWP46Ec/yr59+wra73333cePf/xjvvrVrzJ69Gg++MEPBtL9oQ99iEWLFjFhwoRA2/UmzhhULF3Fy1vejq38NA9uMqdm05+Cmm4bdlUdBts9mKbfuqbPIjIUmA10qmpder6q1onIk3jNm68B7um12a3AKOB7qhrZZe+tv95A9bbWPvNSqVSkfRL93dRDuOWDxwTebsuWLfzmN7/hQx/6UM+8Rx55hBdffJGysjL+8pe/cOWVV/LrX/+aU089lYceeoh//Md/5LXXXqO9vZ1FixZx/fXXc+211/Loo4/ysY99jBtvvDHvfq+66irWrFnD9OnT+9SqPvCBD/Dcc89l3OY973kPjz32WGCP2YijA8JSP/y/aGY3v9w8mPlHjuYH/3pu3vVLrS+9v2yNCKqqqjjuuONKKSlybPdgmn4jko2ILAIW+f9O8adnishK/+8dqnqD//c04FVgM1DRr6ir8bqr+ZaIzPfXOx3vHZyNwE3Rq+9L0p3fLVq0iCFDhjB27FgWLlzIsmXLepZ94Qtf6Ol2/L777uNTn/oUp59+OgCf+MQnuOOOO3jhhRdQVTo7O7n++usRET7ykY/wjW98oyhdUSaTfET9nlApTuQXHj2W/77iPT3/d3R08M3hwwvePttJP27tFUtXceOCOXzmvL53r486KtvdbHuw3YNp+o1INngdcH6i37wj/A94ieUG8uDXbk7lnY4434/XEefdxNARZ6Yax759+xJ9KfKRRx7h7//+7zMuO/zww3v+3rx5Mz/60Y+45553KoD79+9n27ZtdHZ2Mm3atD69JMycOTM+0RETVQ8IcZ+oczUtbmho6LndGdU+4vLzlSdr+cqTtX32FZX+JLHdg2n6jUg2qvol4EsFrlsPZH29WlW3AFdEoSsMpvW02pveyePwww/npptu4qabDqzsPf300zQ2NqKqPds0NDQUPOpfpq58LrzwQp599tmM659zzjk8/vjjBZVdCFHEIK4Tc6HvrkyePDnWfcfhr2Lpqp59xKG/1NjuwTT9A6mBgBHka+llCkuWLOG73/0uL774Yk/vAatWreLtt9/mtNNOY8iQIXzrW9+is7OThx9+mJdeeqnP9iLC6tWrM5Y9efJkXn/99T7zHn/8cdra2jJ+eieadC8A3d3dPb0KBP1Oi4nBlpa9sZyI61csDPSS5K5duyLX0Jugegol3WItbv2lwHYPpul3ySZikn5mUyinnnoq3//+97n22msZN24cRx55JCtXrgS8Zx4PP/wwK1euZPz48fz85z/v09Bgy5YtjBkzJuvDx09+8pNUV1dTXl7OokWLAum67bbbKCsrY8WKFTz00EOUlZVlbL6di7AxqFi6inO+ErbD8cyEPamX6lZsXElnyY/XR15mqbG9j0DT9BvXXY0JFNNdTWdnp9G30gohn4eHHnqIDRs2sHx5zl5/EqOzs5NNmzYF6q4m6tpMsSfw5uZmJk2aFJGawonqezh6bIrXdg+yusubpGIQFUnoz9VdjRHPbAYSqVQq/0qGk8/DZZcFHkC1pASNQZSJJqqTa6HvNUVNWn+x38khw+gpx9aEk1QMosI0/Xbc87GIIUPsz9+2ewii38REA1BeXh5ZWWGoX7GQE6aNDr391j3vNBIpxfg5cZB0DIrFNP12n1UMpLOz08pRLntju4dsQ3b35pibV7EnorYccVy5b9++PfGRRh/9N+9l0jCJ4l3lStPevq0SbavlmBCDYjBNv6vZRMywYcOSllA0tnvIp79iqdmJBmDGjBmxlBuGMI0IXmzO/HaCTTUck2IQBtP0u2QTMabdJw2D7R5y6Y/qZBdXK640GzfG2l9sKIL4PWJM9oZHtiQcE2MQBNP0u2QTMSNHjkxaQtHY7iGb/igTTdyY1KdVb+pXLKSQem9VS+5Tiw0Jx9QYFIpp+l2yiZikhxiIAts9ZNIfxckt7tpMb0zrHr43Gwv4Ho4pz98i0PSEY3IMCsE0/S7ZRIwpQwwUg+0e+uuPKtGUEtO6h89Eru9kw67CTi0mJxwbYpAL0/S7ZBMxttcKwH4PvfXbmGjAvKvSbGT7boIM/mZqwrElBtkwTb9LNhGTr1bQ2R3fS5+9h4VOf6699trA5RRas1m9ejXTp08PrPF3v/tdYE3gjZVz9NFHM2jQoJ6udTKR1m9rogHzrkpzUb9iIV/+x769NQQd/M3EhGNTDDJhmn6XbCIm08Bdu9s7uWPVqxx7yxMcddPjHHvLE9yx6lVa9+V/HyQo6WGh059vf/vbgcuIY/CxKDjhhBP4zne+w8knn5xzvb1797K1pb2ofX3mnIpE3wlZv96uvsU+ftYRfb6vi2cFH9batIRjWwz6Y5p+l2wipqysrM//u9s7+cdvP8fKNW/Q1tGFAm0dXaxc8wYfvOc5drdHn3Ay8ZnPfIYPf/jDPf/feOONzJ8/H1WlpaWFD3zgAxx66KGMGzeOSy65hK1bt/asu3PnTq644gqmTp3KuHHjWLRoEXv27OHCCy9k27ZtPbWobdu25dSwePFiGhoa+OAHP8jo0aP56le/GsjDNddcw/z58/N2MFi3s7jxbOpXLOTGhcFHR42SY45Jdv9hSSecRzeH7wzVFGyNQRrT9LtkEzH93/G49+lNNO1qZ3933/cO9ncrTbvauffpTSXRdeedd1JVVcXKlSt59tln+cEPfsCPfvQjRIRUKsUVV1zB5s2baWhoYNiwYX1uvy1evJi9e/eyYcMGmpub+dznPseoUaN4/PHHmTp1ak8taurUqTk1PPjgg8yYMaOn9vUf//EfgNetRrbPihUrAvl8Zeuunn65wmDKG+6bNpXmuIiD+hULOW9q+A5+TUk4NscAzNPvuquJmP5vr//kpYYDEk2a/d3KT15qYNnCwnsnzkd6WOg0X/va11iyZAkjR47kwQcf5MILL2TMmDHcc889Pc9bJkyY0KfWc9NNN/WM9tnU1MTjjz/OW2+9xbhx4wA499xzI9ML0Y278cpWr5ywvQOYkmiAwM/CTOPrl8/jt7c9E3p7E7q2sT0Gpul3NZuI6T1wV2d3ij37c5/59nR00RVho4FHHnmEXbt29XyWLFnSs+z000/niCOOQFW55JJLeubv3buXT33qU8ycOZNDDjmEefPmsWvXLrq7u9myZQvjx4/vSTSmkk40ACNCHNVJn9j6s2PHjqQlFMWOHTuK/k6TruEMhBiYhEs2EdN74K6hgwcxaljuyuOo4UMYMrg0Ybj33nvp6Ohg6tSpfZ6X3Hnnnbz22mu8+OKLtLa28vvf/x4AVeXwww9n586dGWsfmYZ/zkembXq3nuv/ueOOO/KW2TvRAHQGvINjWqIB7zuxmbR+mxPOQImBKbhkEzH9B6O79N0zGDY480l52GDh0neXprO8jRs3cvPNN/PQQw/x4IMP8tWvfpV169YB8Pbbb1NWVkZ5eTk7d+7sMzLmYYcdxoUXXsjVV19NS0sLnZ2d/PGPfwS84Z/feustdu/e3bP+6tWrcyahTENGZxsuuq2tjWXLlvWst3//fvbt24eq0tnZyb59+1jXsPOAfQQ5qE1MNFBYz9Um01u/rQlnIMXABFyyiZlr5x/JYeVlByScYYOFw8rLuHb+kZHuL93SK/256KKL6Orq4rLLLuPGG2/khBNOYM6cOdxxxx0sXryYjo4Orr/+etrb25k4cSJnnHEGCxYs6FPmgw8+yNChQ5k7dy6TJk3irrvuAmDu3Ll87GMf44gjjqC8vJxt27axZcsWzjrrrKz6vvCFL3DbbbdRXl7O17/+9UDeFixYQFlZGWvWrOGqq66irKyMl19cc8B6+wu8K2lqogH7B+Hrr9/GhDPQYpA0bljoDBQzLHRXV9cBg3ftbu/k3qc38ZOXGtjT0cWo4UO49N0zuHb+kRwywrwhpDN5KJR//dd/5eKLL+b888+PWNWB9L991pvtDa+z5FdNWZebnGgAWlpajH9Olots+otJGsOB10oYt4Eagzhxw0KXkEwn6rFlQ1m28F0sW/guOrtTDC3RM5qwFJNs7r///ojVZCZXohmW5+s1PdGA926TzSe6bPrrVywMnXA6ihUVkIEag6Qw+6xnIfkG7jI90YD5g6flSjQAWVqaA3YkGiDvO0umk0t/MTEo5e20gRyDJDDmzCci00XkARHZJiIdIlIvIneJSKDULCLvEZFH/e33iUiDiPxGRC6IS3tvOjpKff0VPSZ7yJdoAMZkuTNpS6IBeOONN5KWUBT59NuQcAZ6DEqNEclGRGYDLwNXAC8B3wReBz4L/ElEJhRYzmeAZ4H5/vSbwDPAucDjInJT9Or7kq8rFRsw1UMhiQZgd4beamxKNOA1vrCZQvSbnnAOhhiUEiOSDfAdYBJwnaouUtWlqnoeXrI4Grg9XwEiMhRYDuwDTlHVxar6BVVdDJyKd8v3JhEZXozQfC08TO3EMggmeig00QCUD9M+TdBtSzRAT7N0WylUv8kJ52CJQalIvDWaX6vZBNQDs1U11WvZGKAJEGCSqmYdaEVEJgN/A15R1RMyLH8FOA6YqKpv5dKUrTVaQ0MDIsLkyZMZOnRoqJcaHeEoJNmoKqS6aWt5i5fq3+LuF1qsTDQHI8UkDhdjczC9Ndo8f/pk70QDoKpvi8jzwALgDOD3OcppBt4EjhKROapam14gIkcBc4B1+RJNLqZPn86OHTvYvHlzn25perN//37jH7DnwzQPhQ4XoMDezhS7drdy74vtnDR9TLzCYuTll182bjySIATVX0wrtbg42GIQNybUbL4G3ADcoKp3Zlj+beAa4GpV/e88ZV0MPIR3y+yXwDZgGnAR8Brwz6qatyvUbDUbR+lxV7wHD637Ojn+S0+G2tbF2gxy1WxMeGYz1p/uzrI8Pb88X0Gq+gvgPGAX8HFgKbAY2AP8EK/RQUZE5CoRWSsia5uamtixYwdNTU00NjbS0tJCXV0d7e3tVFdXk0qlqKysBN4ZerWyspJUKsXTTz9Ne3s7dXV1tLS00NjYSLq8+vp62traqKmpoaurq2dwo3QZ6WlVVRUdHR3U1tbS2tpKQ0MDzc3NNDc309DQQGtrK7W1tXR0dFBVVZWxjPXr19PV1UVNTQ1tbW3U19cX7Ompp57q46m6ujoRT1+8/1fAO0MMXzyrm8HiXRyVDVbOnJRi9hjl2HEpTpqQYsZo5b1TUlw+p4vffHxWzjgl5anQOD3//POBjz2TPK1duzbwsXfIiKF83I/1pbO7GYTygRndlA9TThifApRDhirHlHtxP3NSikNHKBdMT3HUssci91RZWRnJ7ympOK1evTq2c0Q2T7kwoWZzH7AEWKKqB7wRKCK3A8uAZaq6PE9ZlwHfBx4G/gvYDMwEvghcCvxCVS/JXoJHMTWbVCrVpzNOGzHBQzE1mtfvuDBx/cViQgyKoRj9ptRmD+YYhMX0mk265jI2y/L0/F25CvGfyzwAbAAWq2qNqrarag1e7eZl4GIReV+xgnNRU1MTZ/ElwWYP9SsWWq0/je0eitFfTMKYE+Fzn4M5BnFgQrJ5zZ8elWX5HH+6MU85C4ChwDMZGhqkgD/6/8b6xGzWrFlxFl8SkvYQ9so2fZJKWn8U2O6hWP1hE06U/Rwf7DGIGhOSzR/86QIR6aPHb/p8NrAXeCFPOen3Zw7Nsjw9v7gB6vOwbdu2OIsvCUl6CJto7rrk+J6/XQySJwr9YRNOVK3aXAyiJfFko6p1wJNABV6rs97cCowCHuz9jo2IzBWR/q/HPutPPyIix/deICInAh/Bax37dGTiMzB+/Pg4iy8JSXko5iSx6OTDe/52MUieqPQ//JkzQ20XRcJxMYiWxJONz9V478l8S0QeEZHlIvI08Dm822f9u5l51f/0oKov4bU4KwP+LCI/E5GviMjPgReBEcDdqrohTiMmvn0fFNs89L8Ctk1/Jmz3EJX+k2eGP2EWm3BcDKLFiGTj125OBVYCpwOfB2YDdwNnBHgR85N4/av9CTjfL+cfgOeAj6nq56JVfiA2t15Jk4SHYp/T9MbFIHmi1J/UOzQuBtFijBpV3aKqV6jqYao6TFVnqur1qtqSYV1R1QP6ilGPlar6PlUdp6pDVHW8qs5X1Z+VwsfQoeYNhhaUUnuIMtGAi4EJRK0/iec3LgbRYkyyGSi0tbUlLaFoSukh7MngzIpsLeVdDEwgDv2lTjguBtHikk3ETJw4MWkJRWODh59++j1Zl9mgPx+2e4hL/9Qx4bpzPP223wbexsUgWlyyiZitW7cmLaFoSuUh6ttnaVwMkicu/WtuOj/UdtvbugNv42IQLYl3V2MixXRX09XVxZAhJnSmHZ5SeIgr0YCLgQnErT/O4yeNi0FwTO+uZkCxYUOsLatLQtwevvdMbf6VMnDjgjn5V8LFwATi1l+K5zcuBtHiajYZcEMMxEsprkodA5/3Ln+KhkxjgOfBHUfx4Wo2JSTd9bbNxOmhFInGxSB5SqH/j1/4h1DbFXoMuhhEi6vZZMDVbOLh2C+uoi1ET4nuStSRC1dTNgdXsykhpl1NhCEuD2ESzcSRB7y7mxcXg+Qppf64koaLQbS4mk0GXM0metzVpyNO3PFlBq5mU0LSQ7DaTNQeSn0icDFInlLrj6N1motBtLhkEzFHHZVtDDh7MMFDMVecJugvFts9JKE/7DHzlcczNxF2MYgWl2wipqGhIWkJRROlh6gGsgqCi0HyJKX/+KmjA2/z38/UZ5zvYhAtLtlEzOTJk5OWUDRReUjqPrqLQfIkpf9X150bartMx6qLQbS4ZBMxu3btSlpC0STpIYoHti4GyWPjMfQ/a17v87+LQbS4ZBMxI0aMSFpC0UThIYnbZ2lcDJInaf1hEs5//qrP4L+JeygW0/S7ZOOInKNdM1SHpSR5kTTQcckmYvbt25e0hKIp1kNHiG2iTDQuBsljgv5ijykTPBSDafpdsomY8vLypCUUTTEeTLgyPNhjYAKm6A+TcNLHsCkewmKafpdsImb79u1JSyiasB4+fO8fQ20X9e2zgzkGpmC7/oqlq6z3YJp+l2wiZsaMGUlLKJqwHl7e8nbgbeJ4TnMwx8AUTNIf9hgzyUMYTNPvkk3EbNy4MWkJRRPGgwm3z9IcrDEwCdP0h0k4Sx94IgYlpcO0GLiOODPgOuIMTphk41qfOUqJ66wzflxHnCXEtG69wxDUg2mJ5mCMgWmYqD/oMbf4yO6YlJQG02JgTLIRkeki8oCIbBORDhGpF5G7RGRciLJOFpGfiMhWv6ztIvKMiHw8Du29OeWUU+LeRewE8WDS7bM0B1sMTMRU/d+77KSC131w02DAzGO8EEyLgRHJRkRmAy8DVwAvAd8EXgc+C/xJRCYEKOta4M/AAuD3wJ3AL4HBwPujVX4gpl1NhCFuD3HflnAxSB5T9Z9/7NSC1+1ds/k/v6iMQ06smBYDI57ZiMgTeMnhOlW9p9f8bwCfA76nqp8uoJwFwG+Bp4CPqOrb/ZYPVdW840W6ZzaFYdrtM4ejUNyxGw9GP7PxazULgHrg3n6LbwH2AItFZFQBxX0NaAcu7Z9oAApJNMWyfv36uHcRO4V4WHDn0yVQEo6DJQYmY7r+qWOG5F3n4ll9n9nYdjvNtBgknmyAef70SVVN9V7gJ4zngZHAGbkKEZFjgeOBJ4GdIjJPRG4Qkc+LyHwRKYnXY445phS7iZVCPGx8sz1wuaW6MjxYYmAyputfc9P5edd5dLMJp8fwmBYDE77No/1ptkbhtf4037Bzp/nTZmA18DReTefrwO+AdSJyZHiZhbFp06a4dxE7+TyEucK7ft4RYeUE5mCIgenYoD/fxc9Zkw98xGBT7ca0GJiQbMb6091Zlqfnl+cpZ5I//SRQASz0yz4KeAg4DlglIsMybSwiV4nIWhFZ29TUxI4dO2hqaqKxsZGWlhbq6upob2+nurqaVCpFZaX3wDD9EK6yspJUKsW+fftob2+nrq6OlpYWGhsbSZdXX19PW1sbNTU1dHV19VRz02Wkp1VVVXR0dFBbW0traysNDQ00NzfT3NxMQ0MDra2t1NbW0tHR0TPOeP8y1q9fT1dXFzU1NbS1tVFfX1+wp927d/fxVF1d3cfTSRNSHDsuxewxypmTUpQNVkAZItpz6yH9cDU9nT+1q2SeJk6cmDdO/T2ZFqfhw4cHPvZM8jRlypRQx16pPZ12aIrDRirnTU0xaohy0UzveB0kyp+2S8/xe/GsboYNUi6YnmL+8seM9pSOU3d3d2zniGyecpF4AwERuQ9YAixR1fszLL8dWAYsU9XlOcr5AnCH/+9ZqvqnXssEr5XbqXjPc36aS1MxDQTq6+upqKgIta0p5PJgw4PVgR4DG7BJvw3HdBiSiIHRDQR4p+YyNsvy9PxdecpJL/9b70QDoF5GfdT/990B9QVi9OjgY6CbRjYPx98S/EcpxYoJwUCOgS3YpD9b4hgi2S/EbbidZloMTEg2r/nTbM9k5vjTfB39pMvZlWV5iz8tK0xWODo7Y2/wFjvZPLSGGKjmjQSuAAdyDGzBdv0AU0cm/1pIMZgWAxOSzR/86YL+LcZEZAxwNrAXeCFPOS/gNZOuyNJM+lh/+kYRWvOSSqXyr2Q4mTzYdKthoMbAJmzTn+lYHT889zam125Mi0HiyUZV6/CaK1cA1/RbfCswCnhQVfekZ4rIXBGZ26+cvcAPgBHAbf5zmvT6xwGXA13A/43exTuMHDkyzuJLgu0ebNcP9nuwUX/5iL7/79yf/ybwKV/+TUxqise0GCSebHyuxmuy/C0ReURElovI03i9B2wEbuq3/qv+pz9fBNYB1+N1c3OniDwEvIiXhG7wk1ts7Ny5M87iS0J/DzbVamBgxsA2bNS/7kt9j9mK0flvo72119xbbabFwIhk4yeAU4GVwOnA54HZwN3AGar6VoHltALn4LVKGw9cC3wAeA44X1Xvjlx8P6ZOLbzvJVPp7WH+134ffPsC3s6Ok4EWAxuxVX/vi6RXdhbWvMXU22mmxcCIZAOgqltU9QpVPUxVh6nqTFW9XlVbMqwrqprxSFDVNlW9SVWPUtXhqlquqgtU9cn4XcAbb8T6SKgk9PZQ99a+wNsX8nZ2nAy0GNiI7foB3jPFrGceQTEtBom/Z2Mixbxnk0qlGDTImBweirQH226fpRlIMbAV2/VXLF3FIJRUgMb7Jhz7vUkiBqa/ZzOgWLduXdISisZ2D7brB/s92K5/xthh/PPsYDWbrzy+ISY14TAtBq5mkwE3xIB9jQIcjqhxv4HguJpNCTFtwKIwhPEwZfTgGJSE42CNgUnYrh/g/108JfA2JjUWMC0GrmaTgYO9ZuOu6BwOD/dbCIar2ZSQdK+otnLk0lVcOrs7/4q9MO3HZXsMwH4PtusHz0OYY9uU2o1pMXDJJmJOPPHEpCUURRfwszq7DwvbYwD2e7BdP7zjYUTu1YzFtBjYfVYxkJqamqQlhCZ9Rfb+GYW3wjGtVgN2xyCN7R5s1w/veKixtHZjWgxcsomYWbNmJS2haJ77m92HxUCIge0ebNcPfT2Euai6+qGXopQTGNNiEPisIiKvF/iJtQ8yU9m2bVvSEkLR+0rs+PGFNRoxsVYD9sagN7Z7sF0/FO/hN399MyIl4TAtBmEuYQfhjYnV/zMOr+fmCmBYyLKtZ/z48UlLKJo33s7/1vQT159TAiXhGAgxsN2D7frhQA+2NRYwLQaBE4KqVqjqrAyfcXgDoP0WqAPeFbVYG9i7d2/SEgLT/wcxZFD+ms3RUw6JS07R2BiD/tjuwXb9YL8H0/RHWvtQ1U3Ah4BpwC1Rlm0LtvUHddKtB155Ne7JXbMx9fZZGttikAnbPdiuHzJ7sKl2Y1oMIlejqvuAp4CPRV22DQwdOjRpCYFoaT9w3v5U4Z0PmohtMciE7R5s1w/ZPVSU5xnC0xBMi0Fcqa8LCN7XwwCgra0taQkFk+2Ka+zQ7LfRTK/VgF0xyIbtHmzXD9k9rF7694HLSqJ2Y1oMIk82IjIRuAjYEnXZNjBx4sSkJRTNoCwVG3N6P8vNQIiB7R5s1w+5Pdhw0WVaDMI0ff7PLJ8vi8gP8YZxPgz4fuRqLWDr1q1JSyiIXFdap0/KXLOps+AHBvbEIBe2e7BdP0TvodS1G9NiELgjThHJ93p5K3C3qlrbQKCYjji7uroYMiTZYZELIdeBP2yQHvDc5q5LjmfRyYfHLSsSbIlBLmz3YLt+KMxD0AQye8IIfv9/5hcjq2CSiEHUHXHOy/I5FzgemGhzoimWDRvMGkApE/l+IP8088DrCVsSDdgRg3zY7sF2/RCPhzDDrIfFtBi4IQYyMNCHGAh6NWbD/WmHIyncMATv4IYYKCGmDVjUn0J+GIuPDDbEgGmYHoNCsN2D7fqhcA+mNoQ2LQauZpOBgVqz+cDdq/lr055A2wzUKzCHI0pc7cbD1WxKiGlXE70pNNG4mk3y2O7Bdv0QzMMtC+fGqCQcpsXA1WwyMBBrNu7Ky+GIF/cbczWbklJVVZW0hKK5aGY3Y8zq6SIQAyEGtnuwXT8E92Ba4jAtBsYkGxGZLiIPiMg2EekQkXoRuUtExhVR5ntFpFtEVERui1JvNo466qhS7CYQQa+4nmwcRNV/mfXDCYKJMQiK7R5s1w+l8RDni56mxcCIZCMis4GXgSuAl4BvAq8DnwX+JCITQpQ5BvgRUNJ+thsaGkq5u1i4/e8nJy2hKAZCDGz3YLt+COfBpNqNaTEwItkA3wEmAdep6iJVXaqq5+ElnaOB20OUeTcwFlgencz8TJ5s1ok6zJXTeScfHYOS0mFaDMJguwfb9UPpPMRVuzEtBoknG79WswCoB+7tt/gWYA+wWERGBSjzn/BqSdcBJR0bddeuXaXcXeTUr1hovQfb9YP9HmzXD+E9hKndXP7AmlD7yoVpMUg82eB1dQPwpKr26SdFVd8GngdGAmcUUpiITMLrBPQRVX0oSqGFMGLEiFLvMithr5hM8hAG2/WD/R5s1w+l9bB6Y0vkZZoWAxOSTfqezcYsy2v9aaFPu76P5+vTQUSIyFUislZE1jY1NbFjxw6amppobGykpaWFuro62tvbqa6uJpVKUVlZCbzTlr2yspJUKsXmzZtpb2+nrq6OlpYWGhsbSZdXX19PW1sbNTU1dHV1sX79+j5lpKdVVVV0dHRQW1tLa2srDQ0NNDc309zcTENDA62trdTW1tLR0dHT4qR/GevXr2fYIOWC6SkOHaGcOSnF6CHKIJRjx6WYMVp575QU5cOUD8zoZhDKw5cc5n3htbV9PFVXVxvjqauri5qaGtra2qivr88Yp46OjrxxMt3Tzp07Ax97Jnnq7u7OGyfbPBVy7KU9vX7HhVw623tf7TL/vbV/PqIbUMYMVcYO9X5/M0YrJ01Icey4FO9f8ViknhobGyP1VEiccpH4ezYich+wBFiiqvdnWH47sAxYpqo5n7+IyJXAD4CPqur/+vMuB34I3K6qNxeiqZj3bBoaGpgxY0aobaOkmDb/pngIi+36wX4PtuuH4j2cdfsTbHu7K9A2UTYwSCIGB8V7NiJSAdwF/CKdaJKgvLw8qV0XRe+D3FYPaWzXD/Z7sF0/FO9hzU3nB94mysYCpsXAhGSz25+OzbI8PX9XnnIeANqBqyPQFJrt27cnuXug+APWBA/FYLt+sN+D7fohGg93XXJ8BErCYVoMTEg2r/nTbM9k5vjTbM900pyM13z6Tf8lThURxbuFBnCTP++RotTmIelbB5/9afDbf/2r7kl7KBbb9YP9HmzXD9F4CDMOVFS1G9NiYEKy+YM/XSAiffT4L2aejfdi5gt5yvkfvOc1/T9/9Jev8/9/KhLVWdi4MV9OjJdH1we7msnUK03SHorFdv1gvwfb9UN0HpJ60dO0GCTeQABARJ7Ae9fmOlW9p9f8bwCfA76nqp/uNX8ugKrWFFD25ZSwgUCSnPLl3/DW3mDxNOmNZ4djoHKwdNJpQwOBq4Fm4Fsi8oiILBeRp/ESzUbgpn7rv+p/jCPJbr2DJpq5k0ZmnG9a1+RBsV0/2O/Bdv0QrYckEodpMTCiZgMgIocDXwYuACYATcAvgVtVtaXfugqgqlJAuZdzENRsZi1dRdBI2njl5HDYysFQu7GhZoOqblHVK1T1MFUdpqozVfX6/onGX1cKSTT+uiv99QtKNMWS1NVE0ETz8GfOzLrMtCuioNiuH+z3YLt+iN5DqROHaTEwpmZjErbVbA6GKyaHYyAw0H+rVtRsBgrp7iVMJt/Ba4OHXNiuH+z3YLt+iMdDKROHaTFwySZijjnmmJLuL47uyUvtIWps1w/2e7BdP5jjIexv3BT9aVyyiZhNmzYlLSEnhVxZme4hH7brB/s92K4f4vNQqtqNaTFwySZipk+fXrJ9xTXoUik9xIHt+sF+D7brB7M8hPmtm6QfXLKJnB07diQtISuFXlGZ7KEQbNcP9nuwXT/E66EUtRvTYuCSTcSMHj26JPuJq1YDpfMQF7brB/s92K4f4vcwanCw9YP+5k2LgUs2EdPZ2Zm0hIwEuZIy1UOh2K4f7Pdgu36I38OG24PXbr73TG3+lXxMi4FLNhGTSqXyr1QkcdZqoDQe4sR2/WC/B9v1Q2k8ZOsyKhvLHy+8c03TYuCSTcSMHBns4CkFQe8Pm+ghCLbrB/s92K4fSuPht/8+L/A2H/3uswWtZ1oMXLKJmJ07d8Zafty1GojfQ9zYrh/s92C7fiidh8tOC9Zq7MX61oLWMy0GLtlEzNSpU5OW0IcwrV5M8xAU2/WD/R5s1w+l83Dbh08IvM05y/MPy2VaDFyyiZg33ngjtrJLUauBeD2UAtv1g/0ebNcPpfXwypcWBFp/y+79edcxLQYu2UTM3Llzk5bQQ9i2/CZ5CIPt+sF+D7brh9J6OGREpjFzc3PCl3JffJoWA5dsImbdunWxlFuqWg3E56FU2K4f7Pdgu34ovYegF4e79+VebloM3BADGTBxiIGgycambskdDodH0N/5CIGa5eb81t0QAyUkjgGLgh6Aw4vcn2mDLgXFdv1gvwfb9UMyHoJeJO7LUVcwLQauZpMB02o2rlbjcBw8BP29DwJeN+Q372o2JaSysjLS8oIeeGOLrdYQvYdSY7t+sN+D7fohOQ9BLxaz9RNgWgxczSYDxdRsUqkUgwZFl8OTqNVE7aHU2K4f7Pdgu35I1kMUw0cnod/VbEpITU1NZGUFPeBmjB0WyX6j9JAEtusH+z3Yrh+S9RDFRaNpMXDJJmJmzZqV2L7/+IV/iKScJD1Ege36wX4PtusH+zz0vzg1Tb9LNhGzbdu2SMoJWqs5fmp0Y1dE5SEpbNcP9nuwXT8k76HY2k3S+vvjkk3EjB8/PpH9/uq6cyMrKykPUWG7frDfg+36wU4PvS9STdNvTLIRkeki8oCIbBORDhGpF5G7RGRcgduPEpF/EZGfiEiNiOwRkbdFZK2IfF5EonmgkYe9e/cWXUbQWs3CYyYVvc/eROEhSWzXD/Z7sF0/mOGhmNqNCfp7Y0SyEZHZwMvAFcBLwDeB14HPAn8SkQkFFHMO8BBwPvBX4B7gJ8A04OvAH0RkRPTq+5JE65V7F58WaXm2tyKyXT/Y78F2/WCvh/TFqmn6TVHzHWAScJ2qLlLVpap6Hl7SORq4vYAy/gZcBhymqh/xy/gUcBRQCZwFXBOP/HcYOjR4h3q9CVqrufLMGUXtLxPFekga2/WD/R5s1w/meAhbuzFFf5rEk41fq1kA1AP39lt8C7AHWCwio3KVo6rrVPXHqrq/3/y3gTv9f98XheZctLW1xb2LPvznPx0XeZml9hA1tusH+z3Yrh/s9lCxdJVx+hNPNkB6XNQnVbXPy7B+ongeGAmcUcQ+Ov1pVxFlFMTEiRNDbxu0VrN80TGh95WLYjyYgO36wX4PtusHszyEqd2YpB/MSDZH+9ONWZbX+tOjitjHlf70t0WUURBbt26Nexc9fOyMiljKLaWHOLBdP9jvwXb9YL+HG1b+IWkJfTAh2Yz1p7uzLE/PLw9TuIhcC1wArAMeyLHeVX7LtbVNTU3s2LGDpqYmGhsbaWlpoa6ujvb2dqqrq0mlUj39DqV7Vq2srCSVStHV1UV7ezt1dXW0tLTQ2NhIurz6+nra2tqoqamhq6uL9evX95RRsXQVi4/sBuCimd0MH+R1I1Q+TDnt0BRHj/U+px2a4rCRyr0Lp9DR0UFVVVUfHenp+vXr6erqoqamhra2Nurr6wv2lK5+pz1VV1eH8tR7WlVVRUdHB7W1tbS2ttLQ0EBzczPNzc00NDTQ2tpKbW1tJJ6mTZuWN06mexozZkzgY88kTxUVFaGOPZM8HXnkkZH8nqLyVL9iYZ9zxKghylFjvZtB6fND73NESjW2c0Q2T7lIvG80EbkPWAIsUdX7Myy/HVgGLFPV5QHL/hDwv8CbwNmq+noh2xXTN9r69es54YTgY4qb1LNzWA+mYLt+sN+D7frBTA9BzhMXz+rmF28MLmkv8Kb3jZauuYzNsjw9f1eQQkVkEfAzoBl4X6GJplhKkWge/syZgfcRBNN+YEGxXT/Y78F2/WCmhyCJ4xdvDI5RSXBMSDav+dNsz2Tm+NNsz3QOQEQuBn4BbAfOVdXX8mwSGaUYsOjkmfG+GWzaoEtBsV0/2O/Bdv1gv4f0LbdSDimfCxNuo80GNuE1fZ7du0WaiIwBmgABJqnqngLK+xfgR0AjMC9MjaaUg6cFPRBe+dICDhlhVvt5h8NRWky67d4bo2+jqWod8CRQwYEvXd4KjAIe7J1oRGSuiMztX5aIfAL4H6ABeG+pbp31Ju6roVIkGtuv6GzXD/Z7sF0/2O8hXbMBM2o3iddsoKd2swavF4FHgVeB0/HewdkInKWqb/VaXwFUVXrNmwf8Di+BPgBsybCrXap6Vz49parZuFqNw+EIi4m1G6NrNtBTuzkVWImXZD4PzAbuBs7onWhyMJN3/FyJ1/tA/8/1UerORLqZYRyUKtHE6aEU2K4f7Pdgu36w38NFM7v7/J907caImo1pFFOz6ejoYPjw4XnXM/GqJE2hHkzFdv1gvwfb9YMdHnKdR0YNUfZ0SZ95cZ9HjK/ZDCQaGhqSllA0tnuwXT/Y78F2/WC/h9MnHViRSLJ245JNxEyePDnvOmGe1ZSSQjyYjO36wX4PtusHOzzkqqm8slOyLksCl2wiZteuXZGXWepGAXF4KCW26wf7PdiuH+z30JXKPD+p2o1LNhEzYkTu8dlMr9VAfg+mY7t+sN+D7frBHg/Zaje79ruajSMArqmzw+EIwyCyN/5Konbjkk3E7Nu3L+uyoAF+4vpzipUTilwebMB2/WC/B9v1g10eMtVupoxMQEgOXLKJmPLy8sjKOnrKIZGVFYQoPSSB7frBfg+26wf7PXTneaul1LUbl2wiZvv27Rnnm9azcy6yebAF2/WD/R5s1w/2eehfuzlmnFnvULpkEzEzZsyIpJy4e3bORVQeksJ2/WC/B9v1g/0eXmzO30CglLUbl2wiZuPGA0dCCBrQuy45Pio5ocjkwSZs1w/2e7BdP9jpoXftZsG0LG2fE8Ilm4g57rjjii5j0cmHR6AkPFF4SBLb9YP9HmzXD/Z7+OXmwgZPm12i2o1LNhHTv1vyoLWaGxfMyb9SzNjetbrt+sF+D7brB3s9pGs3vYcYyEVhaxWP64gzA1EOMWByh5sOh2NgEvS8M0KgZnnx5x7XEWcJ6X01FDTgnzmnImI14bD1ii6N7frBfg+26we7PdSvWFhwzQZgXwnqHK5mk4GoajauVuNwOJIi6Pln9FD4638Vdw5yNZsSsn79eiB4oBcdPyUOOaFIe7AV2/WD/R5s1w/2e3j0Y9MDrd/WGZMQH5dsIuaYY44Jtd1dl54SsZLwhPVgCrbrB/s92K4f7PcQRv9Jt8bXMs0lm4jZtGlT4FrNvDnJvcCZiU2bNiUtoShs1w/2e7BdP9jvYdOmTYFvzbe0xyQGl2wiZ/r0YFVXgB9+MrmuaTIRxoNJ2K4f7Pdgu36w30NY/Wff8WTESjxcsomYJff9IdD6x08dHZOS8OzYsSNpCUVhu36w34Pt+sF+D2n9QWs3ja3xPLxxySZimtuDDVj0q+vOjUlJeEaPNi8BBsF2/WC/B9v1g/0eitF/4TeDXTQXgks2EVKxdBVlQwpvSj5nYlmMasLT2Rlzs5SYsV0/2O/Bdv1gv4fe+oPWbl7dvjdqOS7ZRM3QAN/oUzecF5+QIkilzOrALyi26wf7PdiuH+z3UKz+f/n+cxEp8XDJJiLSLdDe6ijsNtrUMUPilFMUI0caNsRfQGzXD/Z7sF0/2O+hv/6gtZvn63ZHKcecZCMi00XkARHZJiIdIlIvIneJyLiA5Yz3t6v3y9nml1uSpiUVowu7jbbmpvNjVhKenTt3Ji2hKGzXD/Z7sF0/2O8hCv2f/Wk0fUSCIclGRGYDLwNXAC8B3wReBz4L/ElEJhRYzgTgT/52dX45L/nlviwiR0Svvm9vAa/szF+zGW/mo5oepk6dmrSEorBdP9jvwXb9YL+HTPqD1m4eXR/daKVGJBvgO8Ak4DpVXaSqS1X1PLxkcTRwe4Hl3AEcBXxDVef75SzCSz6T/P3Eynum5L9PWnmL2X2gvfHGG0lLKArb9YP9HmzXD/Z7iEr/956pjaScxDvi9Gs1m4B6YLaqpnotGwM0AQJMUtU9OcoZDTQDKeAwVX2717JBeDWlmf4+Xs+lKUhHnO/7yu+ob+l4RweKkr12M2owbLjd7GSTSqUYNMiU65Dg2K4f7Pdgu36w30Mu/XF1Emx6R5zz/OmTvRMNgJ8wngdGAmfkKecMoAx4vnei8ctJAU/0218k9E40ADPyPLMxPdEArFu3LmkJRWG7frDfg+36wX4Ppuk3Idkc7U+zDfidrsMdVaJyimJzW/avdFicO46Qk08+OWkJRWG7frDfg+36wX4PufQHfXaz6c3iW6aZkGzG+tNsbtLzy+MsR0SuEpG1IrK2qamJHTt20NTURGNjIy0tLdTV1dHe3k51dTWpVIrKykrue/o1LvMHKLp0djeCMnaYcshQ5b1TUswYrZw0IcWx41LMHqM8+eljaGtro6amhq6urp4uzNODNKWnVVVVdHR0UFtbS2trKw0NDTQ3N9Pc3ExDQwOtra3U1tbS0dFBVVVVxjLWr19PV1cXNTU1tLW1UV9fX5AngKeeegqAyspKUqkU1dXVtLe3U1dXR0tLC42NjaS/o/r6euM8vfjiiwd4Spdhi6fnnnsub5xM9vTnP/851LFnkqf0p9jfU1KeVq9enTNOF8/qZtgg5YLpKQ4dofxduXdj6bCyFCdN8M5f752SonyY8uCjzxTkKRcmPLO5D1gCLFHV+zMsvx1YBixT1eU5ylmG15DgdlW9OcPyJcB9wH2q+qlcmoI8syn03ucg4HU3OJrD4TCIQs9fA+WZTbrGMTbL8vT8XSUqpygunZ15KFabEk366sVWbNcP9nuwXT/Y78E0/SYkm9f8abZnKXP8abZnMVGXE4j+/Zv9rM6Er7Q4TjzxxKQlFIXt+sF+D7brB/s9FKK/lEPRm3BmTHcvusBvotyD3/T5bGAv8EKecl4A2oGz/e16lzMIWNBvf5Hw1A3nUTFueM//759x4Hs2pQxoFNTU1CQtoShs1w/2e7BdP9jvoVD9+c5PkZ2/VDXxD16zZAX+rd/8b/jzv9tv/lxgboZyvuevf2e/+df5839biJ5TTjlFw/L9p6r0knufDL29CezduzdpCUVhu35V+z3Yrl/Vfg9h9Z+57LHQ+wTWapbzqgk1G4Cr8V7I/JaIPCIiy0XkaeBzeLe9buq3/qv+pz/L/PX/XUR+75fzCHC3X/41cRlIM29WGT+/+h/i3k2sbNu2LWkJRWG7frDfg+36wX4PYfWvieldQCOSjarWAacCK4HTgc8Ds/GSxBmq+laB5bwFnAl8CzjSL+d04IfAKf5+YmX8+PFx7yJ2bPdgu36w34Pt+sF+D6bpNyLZAKjqFlW9QlUPU9VhqjpTVa9X1ZYM64qqZuwTRlV3qupn/e2H+eVdqapb43cBe/dGP+hQqbHdg+36wX4PtusH+z2Ypt+YZDNQsLkvpTS2e7BdP9jvwXb9YL8H0/SbpWYAMHTo0KQlFI3tHmzXD/Z7sF0/2O/BNP2J9yBgIiLyJrA55OYTgR0RykkC2z3Yrh/s92C7frDfQxL6Z6rqoZkWuGQTMSKyVrN012ALtnuwXT/Y78F2/WC/B9P0u9toDofD4Ygdl2wcDofDETsu2UTPfUkLiADbPdiuH+z3YLt+sN+DUfrdMxuHw+FwxI6r2TgcDocjdlyycTgcDkfsuGTjcDgcjthxySYPIjJdRB4QkW0i0iEi9SJyl4iMC1jOeH+7er+cbX650+PS7u+3aP0islpENMdnRIz6PyIi94jIsyLS6u/voZBlRRLLgPuMRL+vNdv3/7c4tPv7nSAi/yoivxSRTSLSLiK7ReQ5Eflk/zGoCiivpDGIUn9SMfD3/RW/J/stvoedIvIXEblFRCYELKvkvwNwDQRyIiKzgTXAJOBRoAZ4NzAPb2TQswvpkdo/GNbgjSL6NPBnvDF5/glv6IMzVfV1g/WvBs4Fbs2yym2q2hWF5gz7XgecALQBW/G+tx+r6mUBy4nkuwhKhPrrgXLgrgyL21T168XozLHfTwP/DTThDTzYAEwGPoQ31Pr/Ay7WAk4kScQgYv31JBADf9/7gUqgGu+cMQo4A6+3/G14veNvKaCcRH4HgBmDp5n6IeCgbjnKiWRQtwT1r/YOlURiMA9vSG8B3ufrfiip7yJB/fVAfQLf/3nAB4FB/eZPwTtxK/BhU2MQsf5EYuDve0SW+bf7Hr5jagx69pHEF2fDB288HQXeyHCgjsG7Ut0DjMpTzmi8Ya3bgDH9lg3yD2AFjjBRv79+Ysmmn45QJ+sov4sk9PvbJnaiy6Fpme/nHltiEFa/wTE4wffwlOkxcM9ssjPPnz6pqqneC1T1beB5YCReVTYXZwBlwPP+dr3LSeFdafTeX1REpb8HEfmoiCwVkX8XkQtFZHh0cmMl8u8iIYaLyGUiskxEPisi80RkcIJ6Ov1pIbdQTYxBEP1pTIvBB/3pKwWsm2gMhsRR6ADhaH+6McvyWmAB3nOY3xdZDn45URKV/t78rN//zSJyjar+3xD6Skkc30USTAEe7DfvDRG5QlWfKaUQERkCfNz/97cFbGJUDELoT5NoDETkBry7JWPxnte8By/RrChg80Rj4Go22RnrT3dnWZ6eX16icoIS5X4fxbuCmo5XS5sLLPe3/bmIXBBaZWlIKgZR8kNgPt7JbhRwHN6zwArgcRE5ocR6VgDHAr9R1SfyrYx5MQiqH8yIwQ3ALcD1eInmt8ACVX2zgG0TjYFLNo68qOo3VfUxVW1U1X2q+pqqLgM+j3cMLU9Y4oBHVW9V1adVdbuq7lXVv6rqp/Ee7JYBXyqVFhG5Di/2NcDiUu03KsLqNyEGqjpFVQUv4X0IOAL4i4icHPe+i8Ulm+yks/zYLMvT83eVqJyglGK/9+Pd7z5RRMYUUU7cJBWDUvBdf/reUuxMRK4F7sZrgjtPVXcWuKkRMShCfy5KGgMAP+H9Eu+21wTgfwrYLNEYuGSTndf8abZnKXP8abb7n1GXE5TY96uq+4B0o4dRYcspAUnFoBSkb5/E/v2LyPXAPcBf8U7UQV5kTDwGRerPRcli0B9V3YyXOI8RkYl5Vk80Bi7ZZOcP/nRB/7eM/av4s/GaNL+Qp5wXgHbg7P5X/365C/rtLyqi0p8VETkaGIeXcEwePjf27yJB0i2HIn8puDciciPwTWAd3om6OWARicYgAv25KEkMcjDVn3bnWS/RGLhkkwVVrQOexHv4d02/xbfiXcU8qKp70jNFZK6IzO1XThte65VRHHhP91q//Cc04h4EotIvIrNEZHz/8kXkULwHpgA/05h6EAiCiAz1PczuPT/Md5EE2fSLyLtE5ICrZhGpAL7t/xuqC58CdX0R74H6y8B8Vc16YWFiDKLQn2QMROQoETng1peIDBKR2/F6A1ijqi3+fONiAK67mpxk6NrhVeB0vPbqG4GztFfXDiKiAP4DvN7l9O+u5iXgXbzTXc1Z/oFgnH4RuRzvnvRzeFduO4EZwPvx7vGuBf5BVXdFrd/f/yJgkf/vFOB8X8ez/rwdqnqDv24F3gtrm1W1ol85gb4Lk/SLyJfwHmj/EdiMV5OcDSwERgC/AS5S1f0x6P8EsBLvqvkeMrdkqlfVlbk8+MtKHoOo9Cccg+vxGuE852t7C6/LnXPxGgj8DS+JVufy4C9L5HcAuB4E8n2Aw/Gu4JuA/XgH2l3AuAzrKlnetAfG4z2Y3OyX0wQ8AEw3WT9e886VQBXeQd6Jl3CeBf4NGBaz/i+ldWX51Pdat6L/vLDfhUn68U4qP8VrPbXLj8GbwFN474pIgvoVWG1qDKLSn3AMjsWrPa3Du13dhZc0/+z7G99vfaNikP64mo3D4XA4Ysc9s3E4HA5H7Lhk43A4HI7YccnG4XA4HLHjko3D4XA4YsclG4fD4XDEjks2DofD4Ygdl2wcDofDETsu2TgcDocjdlyycTgcDkfsuGTjcDgcjthxycbhcDgcseOSjcNhMCLyiIioP5Rx/2X/5S/7QRLaHI4guI44HQ6D8ccS+gtel/Jnqupf/Pnz8cYmqQFOU9W9yal0OPLjko3DYTgichbwDN4YJSfjDXK1Dm88odNUdUNy6hyOwnC30RwOw1HVNcAX8caI/x7eyK9TgOtconHYgqvZOBwWICIC/BZY4M/6qapemqAkhyMQrmbjcFiAeleFD/eadVdCUhyOULiajcNhASIyB6jEG454LLABeLeq7ktUmMNRIK5m43AYjogMB36O1zDgo8By4Dhc7cZhES7ZOBzm83XgJOCrqvoUcAvwPPApEbk4UWUOR4G422gOh8GIyEV4z2peBN6jql3+/MPxmj8PAU5S1dcTE+lwFIBLNg6HoYjIDLyEMgg4UVXr+y3/J+AR4M94iWh/iSU6HAXjko3D4XA4Ysc9s3E4HA5H7Lhk43A4HI7YccnG4XA4HLHjko3D4XA4YsclG4fD4XDEjks2DofD4Ygdl2wcDofDETsu2TgcDocjdlyycTgcDkfs/H/ozb1JKjMx8QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Compute the relative L2 error norm (generalization error)\n",
    "\n",
    "x_test = torch.linspace(0, pi, 100000).reshape(-1,1)\n",
    "t_test = torch.ones(100000).reshape(-1, 1)\n",
    "test = torch.cat([x_test, t_test],1)\n",
    "\n",
    "\n",
    "u_test_pred = my_network(test).reshape(-1,1)\n",
    "u_test = exact_solution(x_test,t_test).reshape(-1,1)\n",
    "\n",
    "relative_error_test = torch.mean((u_test_pred - u_test)**2)/torch.mean(u_test**2)\n",
    "print(\"Relative Error Test: \", relative_error_test.detach().numpy()*100, \"%\")\n",
    "\n",
    "# final time ploting:\n",
    "# for time t = 1\n",
    "plt.plot(x_test.detach().numpy(),u_test_pred.detach().numpy(), label='Pred, t=1')\n",
    "plt.scatter(x_test.detach().numpy(), u_test.detach().numpy(), s = 65,label='Exact, t=1')\n",
    "plt.grid(linestyle='dotted')\n",
    "plt.legend(loc='upper left', fontsize = 12)\n",
    "plt.xlabel('x', fontsize = 20)\n",
    "plt.ylabel('u', fontsize = 20)\n",
    "plt.xticks(fontsize=20)\n",
    "plt.yticks(fontsize=20)\n",
    "plt.grid(linestyle='dotted')\n",
    "plt.savefig('first_.png', dpi = 300, bbox_inches = \"tight\")\n",
    "plt.savefig('first_.pdf')\n",
    "# show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "1d2e95f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data taking out for RNN "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "60f3a938",
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "\n",
    "u_sol = torch.zeros(100, 256)\n",
    "x_test = torch.linspace(0, pi, 256).reshape(-1, 1)\n",
    "t_test = torch.zeros(256).reshape(-1, 1)\n",
    "dt = 0.01\n",
    "\n",
    "# # Assuming my_network is your defined neural network model\n",
    "# my_network = YourNetworkModel()\n",
    "\n",
    "for i in range(100):\n",
    "      # Increment t_test by dt\n",
    "    test = torch.cat([x_test, t_test], dim=1)\n",
    "  \n",
    "  \n",
    "    output = my_network(test)\n",
    "    t_test += dt\n",
    "    u_sol[i, :] = output.reshape(256)  # Reshape output to match u_sol[i, :]\n",
    "\n",
    "# Rest of your code...\n",
    "u1 = u_sol.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "b3c76216",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAycAAAFNCAYAAAAaWRC8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAw0UlEQVR4nO3de9BddX3v8c+nCQlCUChxAkI0GYzlJjWQog7tlKKeCbYlbbUe6NFiB5qhlVaP2inWDu2h03O0tjr1QOt5FEZ0VKTowbRNhyLQoVXhEMLFPAQ0XFqCNEgMaECCwe/5Y68Hdnb2s6/r8ltrvV8zzzz7svZev732+n2/6/tbl+2IEAAAAABU7SeqbgAAAAAASBQnAAAAABJBcQIAAAAgCRQnAAAAAJJAcQIAAAAgCRQnAAAAAJJAcQIAAAAgCRQnANAQth+y/UPbu7v+LrX9TtvPdT32gO3fqbq9AAD0ojgBgGb55YhY0vV3Yfb4N+Yek/QWSX9he3WF7QQAYD8UJwDQMhFxh6Stko6rui0AAHSjOAGAlrH9M5JeJWlT1W0BAKDbwqobAADI1bW293bd/wNJP5L0OttPSFogaYmkSyV9u/zmAQAwP/acAECz/EpEHNr198ns8Vuy+4dIOkLSCZL+Z3XNBABgfxQnANAyEbFD0pck/XLVbQEAoBvFCQC0jO3DJf2qpNmq2wIAQDeKEwBolr/v+Z2T/5s9/vq5x9S5Utd3Jf1edc0EAGB/joiq2wAAAAAA7DkBAAAAkAaKEwAAAAD7sH2F7cdsb5nn+f9m+27b37T9dds/ncd8KU4AAAAA9Pq0pLUDnn9Q0s9HxKsl/ZmkmTxmyo8wAgAAANhHRNxse8WA57/edfcWSUfnMV/2nAAAAACYxnmS/imPN6rdnpPDDzs0lr/siKqbUUvP+sCqmyBJWhTPVN0E1MizPpB1BgBQibvuue/xiHhp1e0Y5vWHHRpP/GjvWK+596mnZiV1J9iZiBj70Czbv6BOcfKz4762n9oVJ8tfdoRu+OIVVTejER5adGzVTQCGWvHsvVU3AZgIMRZVI35Ob+mrT/v3qtswiid+tFefec2JY73m1K/d+kxErJlmvrZPkvQpSWdGxM5p3msOh3W1GEELAPL10KJjn/8DqsZ6iCLZfrmkL0t6R0R8K6/3rd2eEwAAUsSGIFL00KJjGYzERGx/QdLpkpba3i7pTyQdIEkR8QlJF0s6XNLf2JakvdPuiZEoTgAAAAD0iIhzhjx/vqTz854vh3W1HKMpAAAASAXFCQAAAIAkUJwAAAAASEJhxYntA23/P9t32Z61/T/6TLPY9hdtb7N966BfoQQAAADQbEXuOdkj6YyI+GlJr5G01vbreqY5T9KuiHilpI9J+nCB7QEAoBBcqQsA8lFYcRIdu7O7B2R/0TPZOklXZrevkfQGZ9ciAwAAANAuhZ5zYnuB7TslPSbp+oi4tWeSoyQ9LEkRsVfSk+pcLxkAAABAyxRanETEcxHxGklHSzrV9omTvI/t9bY32d60c9cTeTYRAAAAQCJKuVpXRDwh6SZJa3ueekTSckmyvVDSSyTt7PP6mYhYExFrDj/s0GIbCwAAAKASRV6t66W2D81uv0jSmyT1/uLfBknnZrffKunGiOg9LwUAAABACyws8L2PlHSl7QXqFEFXR8Q/2L5E0qaI2CDpckmftb1N0vcknV1gewAAAAAkrLDiJCLulrS6z+MXd91+RtKvF9UGAAAAAPXBL8QDAAAASALFCQAAAIAkUJwAAAAASALFCQAAAIAkUJwAAAAASALFCQAAAIAkUJwAAAAASALFCQAAAIAkUJwAAAAASEJhvxBflAV7ntbB928eadqnjjm54NYAQDpGjY3jII4CaJMi4ijGU7viZByDVjASLsa1dcdhub7fcct25fp+aIeyEydxFHnLM5YSRzEJCpC0Nbo4GaTfikmiRd4FyKTzIuFCSj+BEkfRD3EUKUk9jmJ/rS1O+uldgUmyzVZmAh1Xv7aRaJuvCUmUONo+qcZS4mg7NSGOth3FyQDdKzgJthlSTaKj6G07SbYZmp5IiaPNQxxFapoeR6tke62kv5a0QNKnIuJDPc+/QtIVkl4q6XuS3h4R26eZJ8XJiOZWfJJr/dQ5kQ7S/blIsPXS1kRKHK0v4ihS09Y4WibbCyRdJulNkrZLus32hoi4p2uyv5T0mYi40vYZkv6XpHdMM1+KkzExClgPTU2k85n7vCTXtJFMO4ij9UAcRYqIo6U6VdK2iHhAkmxfJWmdpO7i5HhJ781u3yTp2mlnSnEyBUYB09O2ZNqLUcA0kUznd/D9m4mhiSGOEkdTRBytxFGSHu66v13Sa3umuUvSr6lz6NevSjrE9uERsXPSmVKc5IDkWr22J9N+GAWsHsl0NAz0pIE4uj/iaPWIo6NZ/OIDtfKME8Z70dduXWp7U9cjMxExM+as3y/pUtvvlHSzpEckPTfme+yD4iQnJNdqkEyHI7lWg4Q6PuJoNYijwxFHy0cMLcXjEbFmwPOPSFredf/o7LHnRcR31NlzIttLJL0lIp6YplEUJzkjuZaDZDq+rTsOI7GWgIQ6PeJoOYij4yOOloM4mozbJK2yvVKdouRsSb/RPYHtpZK+FxE/lvQBda7cNZWfmPYN0B8dqzgk1Mlt3XEYy69A9Pt8sTyLQxyYHHG0WPT7dETEXkkXSrpO0lZJV0fErO1LbJ+VTXa6pPtsf0vSMkl/Pu182XNSIM5FyR8JIR+M/uWLZFoc4mi+iKH5IY7miziapojYKGljz2MXd92+RtI1ec6TPScFo7Plg5Gq/LE880EfLx7LOB/0+fyxTPNBH0c3ipMSHHz/ZjreFAj+xaHomw79ujws6+nQz4tDHJ0OfRu9KE5KRAccHwG/HCzn8dGfy8dAz2To3+VgOY+P/ox+KE5KRkccHYG+XCzv0dGPq8XyHx39ulws79HRjzEfipMK0CGHI8BXg+U+HP0XdUF/rgbLfTjiKAahOEFyCOzVYvnPj4SaDr6LwejH1WL5z4++i2EoTipC5+yPgA5gVMTR/oijaeB7ACZDcVIhEuu+COTp4LvYH/01TXwv+6LvImX0V4yC4qRidFSkio2cF9BP08b3g1QRR19AP8WoKE6QBAJ4mvheSKioD/prmvheiKMYD8VJAtreaQncAKZFHCWOAmgGipNEtD2xIl1t3uihX9YL3xdSRRwFRkdxgkq1OWDXSRu/JxIq6qKN/bOO+J6A0RRWnNhebvsm2/fYnrX97j7TnG77Sdt3Zn8XF9WeOmjbxhCBGkDeiKNAOtrWH5GPhQW+915J74uIzbYPkXS77esj4p6e6f41In6pwHYAyMHWHYfpuGW7qm5GKUioAIpAHAWGK2zPSUQ8GhGbs9s/kLRV0lFFza8p2tKZGe0DUBTiKADUV5F7Tp5ne4Wk1ZJu7fP0623fJek7kt4fEbOD3uu5p5/W7tvv6PvcklNWT9lSAIO0YdSvqRu2xE0gDW2Io8A0Ci9ObC+R9CVJ74mI7/c8vVnSKyJit+03S7pW0qo+77Fe0npJOvolB887r37Jt46J9+D7N+upY06uuhmFadpo35Z79wx8/sRjF5fUErTdfAXIpK+pY/ycQxyth2HxsxfxtB7qOsgzSQxF/gotTmwfoE5h8rmI+HLv893FSkRstP03tpdGxOM9081ImpGk17xsaYzThu4Vrc6JFtUbN4kOe11dkyyjfukoOpE2ZcAH6Zg0jg56fR1jKXE0HRQk6SmsOLFtSZdL2hoRH51nmiMk7YiIsH2qOufA7CyqTXUqVJo66le30b5pE+k4713HBNs0dRntqzKZzs079RjaZHWKo0XG0H7zII5Wrw5xlIIkbUXuOTlN0jskfdP2ndljfyTp5ZIUEZ+Q9FZJv2N7r6QfSjo7IsbaMzIpEizmU0YyHTTfOiRXRv3Kl1oyrUMMbeogTx1UHUel9GMpcbR8qcVR9FdYcRIR/ybJQ6a5VNKlRbVhFHVIsChHVcm0V52Sa5OkOtqXejIlhqJbKnFUqteAT1MQR5EHfiE+k+KKm2onn1SqhyJsuXdPUgm1W8ptQ/FSjEvz2X37HUm2lzhajpRjVcptS/X7bIpU4xIGozjpwkrcLiknrF4ptpOkWpw6x6K6thuTqVscrUtb8/bQomOrbkLpiEX1RXHSByt0/lLbkK1jgmpzYi1aSqPrTYg/dS6uUkYczUdd25064ijyQnEyj1SSa0qdvQmasIGfUvtT21iqs1RiTp5S+TzE0Xw1JY6m8hmIo/lpYhytmu21tu+zvc32RfNM8zbb99ietf35aedJcTIEK3lzpJKI8pBSYsX0mhxn2FholqbFnaZ9njYjzuTP9gJJl0k6U9Lxks6xfXzPNKskfUDSaRFxgqT3TDtfipMRsMJPJ4VRoaYmoKZ+rjJVParelvjSls9ZFOJocZr6ucpEHG2sUyVti4gHIuJZSVdJWtczzW9LuiwidklSRDw27UwpTkZU5Ypfdaevu6Ynnqo/XwobTXXVtoRKHK2vquNM0ar+fMTRybUtjpbsKEkPd93fnj3W7VWSXmX7a7Zvsb122pkW+SOMjbP79ju4ln/NVJ1wyrLl3j1cy79m2ppQiaP1QxxFqtoaR/tZcNBBk8TWpbY3dd2fiYiZMd9joaRVkk6XdLSkm22/OiKeGLcxc9hzMiY6Qn20JaHOadvnzUNVo+nEEdRF2+JK2z5vnRFHc/F4RKzp+ustTB6RtLzr/tHZY922S9oQET+KiAclfUudYmViFCcTqKJD1PWQhKp2Vbc1wbT1c9cJCbW6ZUAcHU9b40lVn7uuh3ZV0a+Io6W5TdIq2yttL5J0tqQNPdNcq85eE9leqs5hXg9MM1OKkwnRMdLV1oQ6p4rPX9ekWjbixgtYFmkjjrb786eM2FGeiNgr6UJJ10naKunqiJi1fYnts7LJrpO00/Y9km6S9AcRsXOa+XLOCRqFhALUB+efpIk42sE5KIAUERslbex57OKu2yHpvdlfLthzMgWq97SQUF/Ashiu7EMRiBf9sVzSQuzYF8sjLcSLdqA4mVKZHaVux0tzqE+1SKrpIKFiUsTRdqnb913mdglxtD0oTtAIbIj3V+ZyqVtSRToY5EkDcbQ/lgtQLoqTHFDNV4vEgZQRH0bDcqoWcRQpIz60C8UJ0HBsdFSHhAo0A3G0OsTR9qE4yQmdZ19lHeJDwsAkOLQnTcTRfRFH08JyAspBcZKjMhIrG1WYBEm1fGxop4s4ipTV5fy9MvoRcbSdKE5QW2xwj6eM5VWXpIp0sTFSLuLoeFheQPEoTnJGYgVAHJgOy68cbGgjZcSB9qI4QS2RVCfDcgOA6RBHgWJRnBSAah9IV9HHSdP/88FyLBYb2EgZ/b/dKE6Qu6LPOyCpToflB6SP87fSVnQcTf3756ISKBLFSQ0RFIA0MdpXH8RRAEgTxUlB2EgpBqP+6Ut9xA/1QRwtBnE0HyzHYtDvQXECtBBJFQAApIjipEBU//ligxopo7+jDoijSBlxFBLFCdBabKSgDthYQcqIo0D+KE6QK843QMo4CRp1QBwF0GYUJwVj1C8fjE4hZfTzYrF880EcLUZRyzXVIrWoQR76OeZQnNQUI8DIQ9uSKtCNOAoA6aE4AQAAAJAEihMkj0MRkDIORUAdEEeLxfIF8lNYcWJ7ue2bbN9je9b2u/tMY9sft73N9t22Ty6qPVVi4wUpI6miDoijQHPRv9FtYYHvvVfS+yJis+1DJN1u+/qIuKdrmjMlrcr+Xivpb7P/89rz/Wf04I2zkqSVZ5xQSMMBAAAAlK+w4iQiHpX0aHb7B7a3SjpKUndxsk7SZyIiJN1i+1DbR2avHWquSOlGwYKqzW56UCesWVl1M1CCOoz29YuT/RA7kYLZTQ/u9xjxFGiXIvecPM/2CkmrJd3a89RRkh7uur89e2yk4qSfB2+cTTLJ7r79Di05ZXXVzShUEVdoSv2Qo36JtN/jJFeUadSCZL7XpBhDpXbE0SKkHEfni6GDpkk1nm65d49OPHZx1c2onZQHeSaJpU1je62kv5a0QNKnIuJDPc9fIOldkp6TtFvS+p6jpMZWeHFie4mkL0l6T0R8f8L3WC9pvSQdsXjR0OlTLVDQHKMk1PmmTzGxtiGptuWysXkkUw6drU5bLsM9bgzt99oUY2netu44TMct21V1M1qHoqTD9gJJl0l6kzo7EG6zvaGn+Ph8RHwim/4sSR+VtHaa+RZ6tS7bB6hTmHwuIr7cZ5JHJC3vun909tg+ImImItZExJpDFx4w0rwfvHGWlQu5m9304FRJde49gLwVEfOIoShCXjGQWIoiEPf2caqkbRHxQEQ8K+kqdU7JeF7PjoeDJcW0My3yal2WdLmkrRHx0Xkm2yDpN7Ordr1O0pOjnm8yKlYy5CXPRJhHkYNqpXQoQpFxrukDPW3Zo5aKvOMecbRcTe4vTY91E5rv9It92H6X7fsl/YWk3592pkUe1nWapHdI+qbtO7PH/kjSyyUp2wW0UdKbJW2T9LSk3yqiIU09zOvg+zfrqWMaefVlSWkdJ11UAmzyyfMcjlCOspJpKnGU807Gk0ocLbKISOkwrzYcIttETShKnlt80CTbhEttb+q6PxMRM+O+SURcJuky278h6Y8lnTvue3Qr8mpd/ybJQ6YJdU6iKVwKiZWkWj9ljMqlUqCQVOun7ISaQhxF/ZS1dyOVWIrRpbAHugmFyRQej4g1A54f6fSLLlep87MgU2nVL8S3fAVEwjg0AeOqKp4RRzGOsmMbsRTjIJ4NdZukVbZX2l4k6Wx1Tsl4nu1VXXd/UdK3p51pq4oTiRURoyOpIlXEMQBA0SJir6QLJV0naaukqyNi1vYl2ZW5JOlC27PZKRzv1ZSHdEkl/c4JMK6qj5OmUMAgKRyKUCUO78IoqoqjHN6FUTDIM5qI2KjOOeLdj13cdfvdec+zdXtOJFZIDFZlYUJRhGFSiV+ptAPzq3KQp+pYVuX8qx5cw3DEr7S1sjipUttHXDEcSRXzIaF2EEeB+qMfYz6tLU5I8uin6tE+oE6Io+gnlTiaSjuQFuJW+lpbnEisoHnauuOwqpvQKCRV9CJeoQ5Si12ptQfAcK0uToBuqSWx1NpTZ03+VeOqUTQBxWLwLz/Eq3pofXHCigogdanHqdTbVxd13whNdUAl1XYB6K/1xQnSU8VJ2akmr1Tb1WacxAlgGC4ukh4GUeqD4gQAEkZCxSTK3jhOfSAl9fYBeAHFicpP/oy8poWkBUyvznGUc5KAcpW9HcQgT71QnACJK7t4ynPEte7H0AMAOijiURaKk5ojWADNxWgfkB/2kgP1QHGSYSOgnUhWQH6Io+1EHAWQJ4oToAZI/gAAjI9Bk/qhOAEAAACQBIoTtBZ7I5AyRvuA/BH3gfRRnHRhYwAApkMcbRc29gHkjeIESeFXdefHRkD1+I0iAKMin1WPwZJ6ojgBAKBB2CgGUGcUJwAAAACSQHECAAAAIAkUJ2glzt9AyjhOejScA1StusbRurYbaAuKEwAAAABJoDgBsB9OqMU02PMDAJgUxUkPkipSxuEIAACgyShOKsKx0gAAoG3Y/sEwFCcAAABoFI6EyYfttbbvs73N9kV9nl9s+4vZ87faXjHtPClOAAAAAOzD9gJJl0k6U9Lxks6xfXzPZOdJ2hURr5T0MUkfnna+Q4uTPo2Q7dOnnTEAAACAZJ0qaVtEPBARz0q6StK6nmnWSboyu32NpDfY9jQzHWXPydW2/9AdL7L9vyX9r2lmCgAAACBpR0l6uOv+9uyxvtNExF5JT0o6fJqZLhxhmteqs4vm65IOkfQ5SadNM1MAAAAA+XjWB+qhRceO+7Kltjd13Z+JiJkcmzWRUYqTH0n6oaQXSTpQ0oMR8eNCWwUAAACgSI9HxJoBzz8iaXnX/aOzx/pNs932QkkvkbRzmkaNcljXbeoUJz8j6efUORnm74a9yPYVth+zvWWe50+3/aTtO7O/i8dqOQAAAICi3CZple2VthdJOlvShp5pNkg6N7v9Vkk3RkRMM9NR9pycFxFzu3welbTO9jtGeN2nJV0q6TMDpvnXiPilEd6rcZacsrrqJgAAAAB9RcRe2xdKuk7SAklXRMSs7UskbYqIDZIul/RZ29skfU+dAmYqQ4uTrsKk+7HPjvC6m/O41jEAAAAwjpVnnMBvneQgIjZK2tjz2MVdt5+R9Ot5zrPq3zl5ve27bP+T7RMqboukzsoMAACA/HHkCIYZ5bCuomyW9IqI2G37zZKulbSq34S210taL0lHLF5UWgOB1JywZmUp8znx2MWlzAcAAKBbZXtOIuL7EbE7u71R0gG2l84z7UxErImINYcuPKDUdgIAxsMeaADApCorTmwfMfcLkrZPzdoy1aXHAABAOcrak5u3urYbaIvCihPbX5D0DUk/ZXu77fNsX2D7gmySt0raYvsuSR+XdPa0lx4DRkVyQsrY8zAajl0HgOYp7JyTiDhnyPOXqnOpYQAAAACo/GpdAAAgR1zQIg18D8BkKE6AmuBQNAAARschsvVEcdKFlbh6jDQhZZzjAABAsShO0FrsiQDyxyBPu9QtjtatvUAbUZwAQILYyAcAtBHFCQAAAIAkUJxk6jpK+dQxJ1fdBJSAQxEAABhfXbfv2oziBK3GRn+xjlu2q+omoERsBLRTXeJoXdoJtB3FSQW44g9SxhXT0sHGPgCgbShOxAYA0sZoHwCgankeRl72IC3befVCcYLWY+MfmF7ZyZ890GlJPY6m3j4AL6A4ARJGQgUjfsVr4oVFODyzWix/YHKtL05I/OmpIqhTBGBUjNijKHW/gESqcTTVdqFcbO/VR+uLEyBVJFTMST2ppt4+oM7qXrQC42p1cUJCRTeKAWB8xFF0Sy2OptYeVIt4VQ+tLk6QH0Z28kVCzVcTzikgqaIuUolfqbQDwHhaW5xUleg5Xj1tbU9mnMSJcVAwAagb4lb6WlmcsGJikKoLlKrnj3QRuzCOKgcbqo5jVc6fQZ7RVDlYSyxNWyuLEyBVVSd0pC+VpFplO9gDXQ9VxTPiKFBvrStOUknsGKzqkacqkhsJtT7avnFMHAVQd8Sx6dj+SdvX2/529v+wPtO8wvZm23fanrV9wSjv3arihBUR46BYQKqIZaiLsuMocRvjIJZO5SJJN0TEKkk3ZPd7PSrp9RHxGkmvlXSR7ZcNe+NWFSdVa/toax2VlehSSKhV763CeKpKqiRzjKtNcRRokXWSrsxuXynpV3oniIhnI2JPdnexRqw7WlGcrDzjhEYm1CZcHrUOikx4J6xZ2diEyuWli1d2XGtiHG26VAYdio5zqcTRVJZ3XaQwaNvUbcQSLIuIR7Pb/ylpWb+JbC+3fbekhyV9OCK+M+yNF+bXxjSxwtXXiccu1pZ79wyfsAQnrFmp2U0P5v6ewLRWnnGCHrxxtpT5pCCFjRlMZi7m5RlLiaPIS1mxtCjP7F2orTv2O+1jmKW2N3Xdn4mImbk7tr8q6Yg+r/tg952ICNvRbwYR8bCkk7LDua61fU1E7BjUqEbvOUklmaIZ8tzLQUKtv5Q2kouMdYwqIm/E0fpq+hEbLYx1j0fEmq6/me4nI+KNEXFin7+vSNph+0hJyv4/NmhG2R6TLZJ+blijGlmckExRpGkSYqqHcXEoQv0VEffaEEdT29hqy+GQ08TCVONoEdqyPqSEbciRbZB0bnb7XElf6Z3A9tG2X5TdPkzSz0q6b9gbN+6wrlRXqJRGWTG93sQ46DCFtiRRpCGPQxNSjaOYTEqHyPbqjo91j6MM8jRL3Q/zKsGHJF1t+zxJ/y7pbZJke42kCyLifEnHSfqr7JAvS/rLiPjmsDeufXFCEk3Hcct2TXK840ApJ9VudUicbffUMSfr4Ps3V92MUnTHxVGTax1iKYM8zUYcbaclp6zW7tvvqLoZfa084wTpa7dW3YwkRcROSW/o8/gmSednt6+XdNK471274mTxiw+sRRIF0HwpJ9U5xEsAQJ008pyT1DDah5RxKAIAAEgFxQkAIHkM8kyHQYhisXynQ/9GN4oTJI+gXz9tusIMSRUAgPxQnNRUape/RD1R+AEAgJRQnBSMUVUAQAoYjChGUcu1TXuggW6FFSe2r7D9mO0t8zxv2x+3vc323bbZFYB5kVSB9ipqkCfVPdBslCJlRfUbBnMxp8g9J5+WtHbA82dKWpX9rZf0twW2BSUhqdYHBV9+SKoAAOSjsOIkIm6W9L0Bk6yT9JnouEXSobaPLKo9VWCDBUhLqiPlQFkYlMgXyxPIX5XnnBwl6eGu+9uzx4C+SAJIGYMRxWC5Au1Bf4dUkxPiba+3vcn2pp1PP1N1c0ZCB0PKiiz0OLQPSBsDPflgOQLFqLI4eUTS8q77R2eP7SciZiJiTUSsOfygA0tpXMrafGgKyQBoDwZ5gHysePbeqpswMvo9qixONkj6zeyqXa+T9GREPFphe4BWoMArDkkVeWEPZLvx/aPNiryU8BckfUPST9nebvs82xfYviCbZKOkByRtk/RJSb9bVFvKxgZKsdi4BjCtNu+Bloij02r78mt7/0GxFhb1xhFxzpDnQ9K7ipo/qnPcsl3auuOwqpuBPtqeUKVOUj34/s2Fvf+SU1Zr9+13FPb+bcEgD9BexNF2q8UJ8XVCQi0HG9kAMB3i6GRYbkCxKE6AligjoXKcdAeDFNNh+SFVFCblIQ60F8VJjsroSBzn+QKSBABMhziaHgZ50HYUJ0ALsAFSPkb9JsMgzwvK2kglPoyG5VQ+4mg7UZzkhA60L5IqgHERR4H6KKvIJy60D8VJDug41aJAGYzlsz+SKrAv4sRgLB+gPBQnQIOVmVA5Tro/CpTRsJyQKgqT6hEf2oXiZEpldpi6HCddBZIHgCYqs+gnjlaPQR6A4mQqVPJpIbHui+WRDmLFYAzypIO4sS+Wx2Bl9ifiaHtQnKAwjABVh4Q6XNkbqSTW/lgu6SF+dLAc0kO8aAeKkwnRQdJEMqlmGVCIYhLEUaSKXJIu4kYabP+k7ettfzv7f9g8073c9j/b3mr7Htsrhr03xckEqugYHIowujYnlTZ/9jogqWISVRT/xJLyMciDmrlI0g0RsUrSDdn9fj4j6SMRcZykUyU9NuyNKU7GxMbFeKoKtm1MrG38zHVEDOlgkCd9bY0pbf3cdUIcTcI6SVdmt6+U9Cu9E9g+XtLCiLhekiJid0Q8PeyNKU7GQGeolzYlmDZ91jxVtbHa9ljS9s9fJ22LLW37vHkgjrbWsoh4NLv9n5KW9ZnmVZKesP1l23fY/ojtBcPeeGGerWwyOkE9nXjsYm25d0/VzShU1cmUQxEms+SU1dp9+x1VN6N0xNL6aUMclaqNpcTRybQ1jvbzw2dikn661PamrvszETEzd8f2VyUd0ed1H+y+ExFhO/pMt1DSz0laLek/JH1R0jslXT6oUew5GUHVyZRDEaZT9cZ7kZr82dqg6thStrZ93jxVvfHa9FjT9M/XZMSVqTweEWu6/ma6n4yIN0bEiX3+viJph+0jJSn73+9cku2S7oyIByJir6RrJQ3dqKU4GYKVfnpVJ1WpmYmniZ+pjdoSY6r+nAzyTO/EYxc3Lu408TNVoer+VXV8aakNks7Nbp8r6St9prlN0qG2X5rdP0PSPcPemOJkAFb2ZmlSEkrlc6RQeE6r6qQqNT/WNP3ztU0q8WdaqXyOJsTRFBBnSvchSW+y/W1Jb8zuy/Ya25+SpIh4TtL7Jd1g+5uSLOmTw96Yc07mkcpKnsKGU9PU+fjpVJIp8jcXc5p0/HQqcbQpjlu2S1t39P0pgdIRR5GiJsbRVEXETklv6PP4Jknnd92/XtJJ47w3e056LDllNQm1AKmNDNVxL0rd2lsnKQ0CNCX+pPQ5Uvp+m4Q4im4p9bOU4g/Gx56TLqzM7VOH0b9Uk2lqBWeT1Hn0jzjaPnMxKuVYShxtnzrH0bZjz4nS3VuS0ihEHlINwqmO/qXaLpQnxbg0SIrtJY6WJ8WYlWKbUK4U4xIGa/WeE1ZYdOtOYFWOANYhkaa8gTSpp445WQffv7nqZuwn9dE/4ih6pbAnpQ5xtImIo8hDK4uTOiTTpo321U3ZyZVEimG641YKCbYOcbSpUjoxfpCyB3zqFkebOMiTutTiKPprTXFCIk1DXZLqnN5kl1eCrVsSRVqqHAWsUyxlkCcdRcRS4miaUt170ou9KelqbHFSpwTai4SarkHJsDfZNjVxNnm0ry5JdU5vnCsiydY5ljZZ3QZ6eg2Lj1vu3dPYGCo1O47WTRlxFOOpbXFCwqyvuifV+TQ5kaIe5ouLw5JtE+Mpgzz1RjxFVZoYD+umdsXJgoMOYsUBKtKG0b667T0ZBTGzmZo60NN0xFFgMC4lnJi2jPa1ITgDqEZb4igANBHFSUJIqEhZmwpK+iLqok39sgna9H0RRzEpihNUpk1Buu74rlAXbBAB6aA/YhIUJ4loawdmoxepamufRP0QR+uB7wkYDcUJgIHanFApUOqlzd9Xm/tpHbT5+2lzv8RkKE4S0PaO2+agDSAfbY+jANAUFCcVI6F2UKCkie+FPor6oL+mie+FOIrxFFqc2F5r+z7b22xf1Of5d9r+ru07s7/zi2wPgNGRUF9AYk0b388L6Ldp4ft4Af0UoyqsOLG9QNJlks6UdLykc2wf32fSL0bEa7K/TxXVnhTRUfdFEE8H38X+6K9p4nvZH/03DXwP+6O/YhRF7jk5VdK2iHggIp6VdJWkdQXOr1booP0RzKvHd4C6II4iVcRRYHJFFidHSXq46/727LFeb7F9t+1rbC8vsD3JIKEORlBHqui7mM+KZ++tugn7II4iVcRRDFP1CfF/L2lFRJwk6XpJV/abyPZ625tsb/ru93eX2sC80SlHQ2KtBst9OPpwGvgehqM/l++4ZbtY7iOg/2KQIouTRyR17wk5OnvseRGxMyL2ZHc/JemUfm8UETMRsSYi1rz0xUsKaWwZ6IzjIcCXh4Q6HvpytVj+o6Nfl4dlPR76MeZTZHFym6RVtlfaXiTpbEkbuiewfWTX3bMkbS2wPZWiE06GYF88lvFk6NPVYLmPj8GH4rF8J0N/Rj+FFScRsVfShZKuU6fouDoiZm1fYvusbLLftz1r+y5Jvy/pnUW1p0p0vukQ9IvDsp0OfbtcLO/p0N+LwXKdDv0avRYW+eYRsVHSxp7HLu66/QFJHyiyDVWj0+VjLvhv3XFYxS1pDhJqPp465mQdfP/mqpvRaMTR/By3bBdxNEfE0XwQR9Gt6hPiG42Emj8SQT5Yjvl66piT6e8FYbnmj/4/PQ6Vyx9xtF5s/6Tt621/O/vfd9TD9odtb8n+/uso701xUgA6WLFICpNj2RWLfp8vlmdxiAWTY7kVi35fGxdJuiEiVkm6Ibu/D9u/KOlkSa+R9FpJ77f94mFvXOhhXW1EpyoPhyeMjmRanrkYwCEKkyOOlodDZkdHHC0PcbQW1kk6Pbt9paR/kfSHPdMcL+nm7Dz0vbbvlrRW0tWD3pg9Jzlhb0k1GP0bjuVTDeLB+Iij1SGWzo9lUx3iQdKWRcSj2e3/lLSszzR3SVpr+yDbSyX9gvb9mZG+2HMyJTpOGhj92x/JtHqM/o2OWJoGYukLiKFpII6O5odP7dHspgfHfdlS25u67s9ExMzcHdtflXREn9d9sPtORITt6J0oIv7Z9s9I+rqk70r6hqTnhjWK4mRCJNI0tT2xkkzTRHKdH7E0TW2OpcTRNBFHC/F4RKyZ78mIeON8z9neYfvIiHg0+93Cx+Z5jz+X9OfZaz4v6VvDGkVxMiYSaT10J5c2JFeSaT10x482J1jiaH20JZYSQ+uDIiUZGySdK+lD2f+v9E5ge4GkQyNip+2TJJ0k6Z+HvTHFyQhIpPXW1ORKMq23NiZYYmm9NS2WEkPrjcGeyn1I0tW2z5P075LeJkm210i6ICLOl3SApH+1LUnfl/T27OT4gShO5kESbabeZFSnBEsibaamJ1hiaTPVMZYSQ5ur6XE0RRGxU9Ib+jy+SdL52e1n1Lli11goTjIk0Hbql6xSSLIk0XbqjUN1TLLE0nZKLZYSQ9urCXG07VpZnJA8McigpJZXsiVxYhT9YlVKiZZYikGGxblp4ikxFKNKPY5if40tTkiaKAIJEVUbFtvyTLrEURSJeIqqDIptFC7Vq11x8tzig0iYADAP4iMATI4YWj1+IR4AAABAEihOAAAAACSB4gRA0h5adGzVTQAAACWhOAEAAACQBIoTAAAAAEmgOAEAAACQBIoTAAAAAEmgOAEAAACQBIoTAAAAAEmgOAEAAACQBIoTAAAAAEmgOAEAAACQBIoTaMWz91bdBAAAAIDiBAAAAEAaKE4AAAAAJIHiBAAAAEASKE4AAAAAJIHiBAAAAEASKE4AJO+hRcdW3QRgKK58CADTozgBAAAAkASKE0hixA8AgCYiv6NuCi1ObK+1fZ/tbbYv6vP8YttfzJ6/1faKItsDoL44tAt1wIYggDaw/eu2Z23/2PaaeaZZbvsm2/dk0757lPcurDixvUDSZZLOlHS8pHNsH98z2XmSdkXEKyV9TNKHi2oPhiOpInUUKKiDFc/e+/wfUCXWQRRoi6Rfk3TzgGn2SnpfRBwv6XWS3tWnFtjPwnza19epkrZFxAOSZPsqSesk3dM1zTpJf5rdvkbSpbYdEVFguzBAv0DGBiFS0r0+kniRukHrKLEVeSMmoiwRsVWSbA+a5lFJj2a3f2B7q6SjtG8tsJ8ii5OjJD3cdX+7pNfON01E7LX9pKTDJT1eYLswJoLdC9iYSMuw74N1Fylj/QTQFtmpG6sl3Tps2iKLk9zYXi9pfXZ3z9JXn7alyvYgSUtFUYv9sV6gH9YL9MN6gX5+quoGjOLJx2ev+4eZ45eO+bIDbW/quj8TETNzd2x/VdIRfV73wYj4yqgzsb1E0pckvScivj9s+iKLk0ckLe+6f3T2WL9pttteKOklknb2vlG2oGYkyfamiOh74g3ai/UC/bBeoB/WC/TDeoF+ejbekxURawt4zzdO+x62D1CnMPlcRHx5lNcUebWu2yStsr3S9iJJZ0va0DPNBknnZrffKulGzjcBAAAA6s2dE1Iul7Q1Ij466usKK04iYq+kCyVdJ2mrpKsjYtb2JbbPyia7XNLhtrdJeq+k/S43DAAAACAdtn/V9nZJr5f0j7avyx5/me2N2WSnSXqHpDNs35n9vXnoe9dtR4Xt9d3HwwES6wX6Y71AP6wX6If1Av2wXpSvdsUJAAAAgGYq9BfiAQAAAGBUyRYnttfavs/2Ntv7nYtie7HtL2bP35pdPxkNN8J68V7b99i+2/YNtl9RRTtRrmHrRdd0b7EdtrkiTwuMsl7YflsWM2Ztf77sNqJ8I+SRl9u+yfYdWS4Zeow86s32FbYfs933pyrc8fFsnbnb9sllt7FNkixObC+QdJmkMyUdL+mcPj93f56kXRHxSkkfk/ThcluJso24XtwhaU1EnCTpGkl/UW4rUbYR1wvZPkTSuzXCD0Ch/kZZL2yvkvQBSadFxAmS3lN2O1GuEePFH6tzEZ/V6lxp9G/KbSUq8GlJgy7Fe6akVdnfekl/W0KbWivJ4kTSqZK2RcQDEfGspKskreuZZp2kK7Pb10h6Q3bJMjTX0PUiIm6KiKezu7eo8/s6aLZR4oUk/Zk6gxjPlNk4VGaU9eK3JV0WEbskKSIeK7mNKN8o60VIenF2+yWSvlNi+1CBiLhZ0vcGTLJO0mei4xZJh9o+spzWtU+qxclRkh7uur89e6zvNNlli5+UdHgprUNVRlkvup0n6Z8KbRFSMHS9yHbBL4+IfyyzYajUKPHiVZJeZftrtm+xnfuPmCE5o6wXfyrp7dllUjdK+r1ymoaEjbv9gSkU+QvxQGVsv13SGkk/X3VbUC3bPyHpo5LeWXFTkJ6F6hymcbo6e1lvtv3qiHiiykahcudI+nRE/JXt10v6rO0TI+LHVTcMaINU95w8Iml51/2js8f6TmN7oTq7XneW0jpUZZT1QrbfKOmDks6KiD0ltQ3VGbZeHCLpREn/YvshSa+TtIGT4htvlHixXdKGiPhRRDwo6VvqFCtorlHWi/MkXS1JEfENSQdKWlpK65CqkbY/kI9Ui5PbJK2yvdL2InVOSNvQM80GSedmt98q6cbgR1uabuh6YXu1pP+jTmHC8ePtMHC9iIgnI2JpRKyIiBXqnIt0VkRsqqa5KMkoeeRadfaayPZSdQ7zeqDENqJ8o6wX/yHpDZJk+zh1ipPvltpKpGaDpN/Mrtr1OklPRsSjVTeqqZI8rCsi9tq+UNJ1khZIuiIiZm1fImlTRGyQdLk6u1q3qXMS09nVtRhlGHG9+IikJZL+Lrs+wn9ExFmVNRqFG3G9QMuMuF5cJ+m/2L5H0nOS/iAi2APfYCOuF++T9Enb/12dk+PfyeBns9n+gjoDFUuzc43+RNIBkhQRn1Dn3KM3S9om6WlJv1VNS9uBX4gHAAAAkIRUD+sCAAAA0DIUJwAAAACSQHECAAAAIAkUJwAAAACSQHECAAAAIAkUJwDQYLYPtf27VbcDAIBRUJwAQLMdKoniBABQCxQnANBsH5J0jO07bX+k6sYAADAIP8IIAA1me4Wkf4iIE6tuCwAAw7DnBAAAAEASKE4AAAAAJIHiBACa7QeSDqm6EQAAjILiBAAaLCJ2Svqa7S2cEA8ASB0nxAMAAABIAntOAAAAACSB4gQAAABAEihOAAAAACSB4gQAAABAEihOAAAAACSB4gQAAABAEihOAAAAACSB4gQAAABAEv4/cUI0anN/ZoEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(0, pi, 256)\n",
    "t = np.linspace(0, 1, 100)\n",
    "\n",
    "X, T = np.meshgrid(x, t)\n",
    "# # Define custom color levels\n",
    "#c_levels = np.linspace(np.min(u), np.max(u), 100)\n",
    "u1 = u1.detach().numpy()\n",
    "# Plot the contour\n",
    "plt.figure(figsize=(15, 5))\n",
    "plt.contourf(T, X, u1.T, cmap='coolwarm')\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('x')\n",
    "plt.title('EB')\n",
    "plt.colorbar()  # Add a colorbar for the contour levels\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "b15f0db6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.io\n",
    "\n",
    "# Example data\n",
    "data = {\n",
    "    'x': x,\n",
    "    't': t,\n",
    "    'u1': u1,\n",
    "}\n",
    "\n",
    "# Save data to .mat file\n",
    "file_path = 'EB.mat'\n",
    "scipy.io.savemat(file_path, data)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9524816a",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pytorch",
   "language": "python",
   "name": "pytorch"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
